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ABSTRACT OF THE DISSERTATION 
Properties of Neutron Star Critical Collapses 

by 

Mew-Bing Wan 
Doctor of Philosophy in Physics 
Washington University in St. Louis, 2009 
Professor Wai-Mo Suen, Chairperson 

Critical phenomena in gravitational collapse opened a new mathematical vista 
into the theory of general relativity and may ultimately entail fundamental physical 
implication in observations. However, at present, the dynamics of critical phenom- 
ena in gravitational collapse scenarios are still largely unknown. My thesis seeks to 
understand the properties of the threshold in the solution space of the Einstein field 
equations between the black hole and neutron star phases, understand the properties 
of the neutron star critical solution and clarify the implication of these results on 
realistic astrophysical scenarios. We develop a new set of neutron star-like initial 
data to establish the universality of the neutron star critical solution and analyze the 
structure of neutron star and neutron star-like critical collapses via the study of the 
phase spaces. We also study the different time scales involved in the neutron star crit- 
ical solution and analyze the properties of the critical index via comparisons between 
neutron star and neutron star-like initial data. Finally, we explore the boundary of 
the attraction basin of the neutron star critical solution and its transition to a known 
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set of non-critical fixed points. 
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Chapter 1 
Introduction 

1.1 Background and motivation 

The field of numerical relativity has been a rich interplay of different disciplines 
throughout its relatively young life. Besides the advances in numerical methods ap- 
plied to the theory of general relativity, eg. spectral methods [S], finite-element 
methods [2], [50] and multigrid methods [6T],[42j, diverse fields such as theoretical 
computational science, eg. adaptive mesh refinement algorithms and algebraic com- 
putation; fluid dynamics, eg. magnetohydrodynamical flows; the theory of partial 
differential equations; the theory of the global structure of spacetimes, eg. the study 
of the structure of singularities and the dynamics of the spacetimes surrouding it; 
theoretical astrophysics, eg. binary pulsars and black holes, rotating stars, and su- 
pernovae core collapses; cosmology, eg. galactic clusters and collisions; and astron- 
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omy, eg. gravitational radiation detection, gamma-ray burst modelling; have con- 
tributed to its development [2E] , [3D] , PS] , [ZSJ , , p3] , pi] , [78] , [77] , [67] , [57] , [7j , [HI] , [35] , 
and vice versa. An example is the application of the Godunov method in dealing 
with hydrodynamic shocks when solving the Einstein field equations on a numerical 
finite-differenced grid. Godunov first suggested this method in 1959 as a means to 
solve partial differential equations. In the application of his method in numerical 
relativistic hydrodynamics, the relativistic fluid is cast as a Riemann problem, which 
consists of a conservation law together with a piecewise constant data having a single 
discontinuity. The use of the Godunov method is now widespread and known to be 
highly effective in all state-of-the-art numerical general relativistic hydrodynamics 
simulations of astrophysical systems and phenomena [31j,|32j. 

In 1993, a phenomenon is found by Choptuik [20] in general relativistic simula- 
tions, analogous to phase transitions in condensed matter systems, thereby earning 
the name critical phenomenon. Critical phenomena are the first new phenomena dis- 
covered via numerical simulation in the theory of general relativity. The phenomena 
observed describes the existence of a threshold between distinct possible end states of 
a gravitational collapse scenario, and new behavior of the gravitating system on this 
threshold reflective of that undergone by condensed matter systems during phase tran- 
sitions, eg. liquid-gas transitions and ferromagnetic phase transitions. The dominant 
features of this behavior include fine-tuning, universality and scale-invariance, where 
a term called order parameter is used to describe a physical quantity that follows 

2 



Chapter 1 



Introduction 



a power-law behavior, and where the phase transitions are categorized according to 
whether the order parameter changes discretely or continuously. The scale-invariant 
feature of these condensed matter systems, eg. the dimensionless quality of the corre- 
lation length which causes the particle clusters to exhibit a fractal structure, and the 
renormalization group transformation performed on the correlation length, is carried 
over analogously to the gravitational collapse behavior on the threshold. As such, the 
discovery of critical phenomena in gravitational collapse became another prominent 
example where a separate field of physics has enriched the field of numerical relativity. 
In addition, it has sinced opened a new mathematical vista into the theory of relativ- 
ity. In the first decade after this discovery, research has been focused on analytically 
establishing the existence and generality of critical phenomena in general relativity. 
Due to the complexity of Einstein's field equations, such analytic analyses of the 
phenomena have all assumed spherical symmetry, and involved simplified matter or 
massless systems. Overview of these works has been covered comprehensively by 
Gundlach in his review [38J. Later works gradually ventured into more relaxed sym- 
metry, ie. axisymmetry, and considered more complicated matter models. However, 
the field of critical gravitational collapse still dwells very much in the theoretical and 
mathematical arena, and has been used only as a tool to understand new and inter- 
esting aspects of the theory of general relativity, eg. the cosmic censorship conjecture 
and the emergence of fractality. 

In 2003, Noble and Choptuik [63] began looking into how critical phenomena could 
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be triggered in static neutron star models in spherical symmetry. Neutron stars refer 
to stars with masses on the order of 1.5 M , radii of about 12 km, and central 
densities 5 to 10 times the nuclear equilibrium density of about 0.16 fm 3 , thereby 
making them some of the densest massive objects found in the astrophysical realm 
[ID] , [H] , [43j . They were first predicted by Baade and Zwicky in 1933 [8] in a study 
on the origins of supernovae, and discovered by Bell and Hewish in 1967 as a source 
of regular radio pulses in the Crab Nebula pflj . Neutron stars are now thought to be 
formed as an aftermath of the gravitational collapse of the core of a massive star of 
more than 8 M at the end of its life [22] ■ They are also believed to be created from 
the accretion-induced collapse of massive white dwarfs. Their interiors, in particular 
their outer cores, consist of neutron superfluids with proton superconductors. These 
interiors lose energy at a rapid rate via neutrino emission. A standard neutrino 
cooling scenario called the Urea process [71] requires the proton to neutron ratio 
to exceed 1/8. Each reaction in this process produces a neutrino and antineutrino 
via alternate beta and inverse-beta decays, losing thermal energy in the process. The 
period whereby this occurs beginning from the explosion is about 3 x 10 5 years old. In 
1974, Hulse and Taylor discovered a binary neutron star system |45j . a system which 
is also termed as a binary pulsar, due to their emission of radio pulses as first observed 
by Bell and Zwicky. Two stars starting from large separations slowly inspiral into 
each other due to loss of angular momentum and energy. When their orbits shrink 
beneath the innermost stable circular orbit, they enter a coalescence phase, where 
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they begin to plunge and merge. Even though this phase is characterized by very 
strong and dynamical gravitational and hydrodynamical processes, the plunging can 
be approximately described as two stars colliding head-on into each other. Due to 
the pervasiveness of such systems in galactic clusters as well as the dynamics of their 
orbital decrease and energy loss, coalescing binary neutron stars are also generally 
considered good candidates as sources of gravitational radiation detectable by both 
existing ground-based detectors such as LIGO, VIRGO, GEO600, and TAMA300 and 
the proposed space-based detector LISA [26J. 

In current state-of-the-art numerical simulations of such processes, full general 
relativistic hydrodynamic is employed, where the standard Tolman-Oppenheimer- 
Volkoff models are solved together with a polytropic equation of state. The study 
of the time scales involved in gravitational collapses in colliding neutron stars was 
first considered by Miller, Suen and Tobias [58] using fully general relativistic hy- 
drodynamics simulations. In this study, the neutron stars are modelled such that 
they infall into each other from infinity. As a result, a dividing line is found to exist 
between neutron star masses producing collapses that occur within a dynamical time 
scale, ie. prompt collapses, and those that occur after neutrino cooling settles in, ie. 
delayed collapses. As the problem is studied using a 3-dimensional code, the quest of 
determining the dividing line between these two scenarios becomes computationally 
expensive. This drove the need to construct a 2-dimensional code capable of high 
resolutions in the finite-differencing scheme employed in the code (refer to Appendix 
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A), which in turn can be employed to resolve the fine structure behavior at the di- 
viding line. In 2007, such an axisymmetric construction was undertaken by Jin and 
Suen jlT] in a study that explores this dividing line in fine detail. They report that 
critical phenomena is observed along this dividing line, characterizing an oscillatory 
threshold between the black hole and neutron star end states. Furthermore, the same 
phenomena is seen to occur when the equation of state of the neutron star system is 
made to vary infinitesimally. Although the same equation of state is used in both this 
study and in the study undertaken by Noble and Choptuik earlier, the former does 
not exhibit scale-invariance at the threshold. With all these considerations, this study 
poses a significant departure from previous works and approaches the realm of real- 
istic astrophysical phenomena. Furthermore, the determination of the gravitational 
radiation signature of the unstable modes of such gravitational collapses of neutron 
star systems may provide insights to gravitational radiation emission data. However, 
at present, the dynamics of critical phenomena in such astrophysical gravitational 
collapse scenarios are still largely unknown. 

In addition, the study of non-rotating neutron stars is still a deviation from real- 
istic systems typically formed in the universe. In particular, coalescences of binary 
neutron star inspirals produce hypermassive neutron stars that support themselves 
against gravitational collapse via differential rotation. These hypermassive neutron 
stars undergo diverse astrophysical mechanisms, eg. angular momentum transport, 
magnetic breaking, as well as gravitational radiation emission, before collapsing into 

6 



Chapter 1 



Introduction 



black hole states. Rotating neutron stars are also formed from rotational core collapse 
of supernovae. With this in mind, Jin and Suen [48J recently incorporated angular 
momentum into the head-on colliding neutron star system. The axis of angular mo- 
mentum is parallel to the axis of collision. In this recent study, they find similar 
critical phenomena in the simultaneously oscillating and rotating merger. 

Given this background, my thesis is motivated by a three-fold objective as follows: 
(1) understand the properties of neutron star critical collapses, (2) understand the 
properties of the threshold in the solution space of the Einstein field equations between 
the black hole and neutron star end states, and (3) clarify the implication of the results 
on realistic astrophysical scenarios. 

1.2 Overview of thesis 

In Chapter 2 of the thesis, the theoretical basis of numerical relativity is presented, 
namely detailed derivations of the 3+1 ADM formalism in solving the system of 
Einsteins field equations in the theory of general relativity. Chapter 3 presents the 
various constructs required in the implementation of the 3+1 formalism in neutron 
star simulations, where we see the coupling of the spacetime with the hydrodynamical 
matter equations adapted to the 3+1 formalism. It further explores in particular how 
conformal decomposition and the BSSN formalism is employed in the solution and 
evolution of the system of Einsteins field equations in the theory of general relativity 
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coupled with hydrodynamically-described matter, the basis of the GRAstro-3D code 
[31], [32]. We also present the main concepts behind the set-up and mechanism of the 
GRAstro-2D code. In Chapter 4, we study the Tolman-Volkoff-Oppeheimer spacetime 
and the ingredients of stellar perturbation theory as well as numerical observations in 
stellar phase transitions in order to enable us to understand the time scales involved 
in neutron star critical collapses. Chapter 5 in turn presents the basic concepts 
of critical gravitational collapse as well as some features of its implementation in 
numerical simulations of neutron stars. In this chapter, the dynamical systems picture 
is brought forward prominently as an important tool in understanding the structure 
of critical solutions. 

Chapter 6 contains the bulk of the analysis and simulation results in line with the 
three-fold objective mentioned above, utilizing all the ingredients presented in former 
chapters. In Section 6.1, curve fitting is employed in order to establish the time scales 
and oscillation frequencies exhibited by the critical solution and their universality 
across several crucial parameters in the critical solution. These results are compared 
with the values obtained from analytic perturbative analysis of equilibrium TOV 
configurations as presented in Chapter 4 of the thesis. The construction of a neutron 
star-like initial data is carried out and evolution of various sets of this new initial 
data is studied and analyzed in Sections 6.2 and 6.3 in order to provide evidence that 
the neutron star critical solution is a semi-attractor. The boundary of its attraction 
basin and the transition to a known set of non-critical fixed points is explored. In 
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Section 6.4, the evolutions of the neutron star initial data sets obtained by Jin and 
Suen [47] are studied in detail in the aspects of its spacetime and matter properties. 
Convergence tests and various parameter searches for the construction of a phase 
diagram are also performed in this section. The power of the radiation emitted at the 
brink of critical collapse is examined using the quadrupole approximation in order to 
shed light on the structure of the semi-attractor, ie. on whether it is a limit cycle or 
a limit point. In Section 6.5, convergence tests are performed for the critical indices 
obtained from different sets of the new initial data in order to determine whether the 
neutron star critical index has a 1-parameter or 2-parameter dependence. The mass 
dependence of the critical index is also studied. 

Chapter 7 summarizes the main conclusions drawn from the analyses and sim- 
ulation results and propose interpretations of them in the framework of astrophysi- 
cal relevance. Finally, the Appendix contains explanations of the finite-differencing 
scheme, of the basic geometric constructs of the pushforward and pullback operators 
used in the previous chapters, and of the derivation of geometric units used in this 
thesis. 
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Chapter 2 

3+1 formalism in numerical 
relativity 

Throughout the past half century, the theory of general relativity has proven to be 
extremely successful in the understanding of phenomena of massive objects interact- 
ing under strongly dynamical gravitational fields. In particular, numerical advances 
within the field have made enormous headway in solving the celebrated two-body 
problem in geometrodynamics. As a result, we are now able to understand to a great 
extent the various physical phenomena observed in binary neutron star coalescences 
and binary black hole coalescences. This includes the understanding of the broader 
context of gravitational collapse of compact objects, the various mechanisms that 
trigger their occurence and their effects on the spacetime structure surrounding these 
systems via the emission of gravitational waves. 
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Figure 2.1: 

Numerical advances in the theory are characterized by the success of fully cou- 
pled general relativistic simulations of the astrophysical systems as mentioned above. 
These simulations employ massive numerical algorithms that solve the Einstein field 
equations coupled with matter. Their theoretical basis is the 3+1 formalism which 
reduces the Einstein field equations to an initial- value problem [5].[85].[S5]. In this 
chapter, we present this formalism as it is adapted to the simulations performed for 
the study in this thesis, ie. the study of critical phenomena in gravitational collapse 
of neutron stars. 
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2.1 Hypersurfaces and foliations 

In the 3+1 formalism, the 4- dimensional spacetime is foliated by 3-dimensional spatial 
hypersurfaces of constant coordinate time. We denote the spacetime by M and the hy- 
persurface by H . The hypersurface H is obtained by an embedding of a 3-dimensional 
manifold in the 4-dimensional spacetime via a 1-to-l mapping that ensures that the 
hypersurface does not intersect itself. The embedding induces a push-forward map- 
ping between vectors on H to vectors on M, and a pull-back mapping between linear 
forms on H and those on M (refer to Appendix B). 

As the hypersurfaces are labeled by constant coordinate time, t, the gradient 1- 
form dt is timelike and can be written in component form as (dt) M = \/^t. We 
denote the future-directed timelike unit vector normal to a hypersurface as n. We 
also define as Eulerian observers a class of observers whose 4- velocity is collinear with 
this vector and thus whose worldlines are orthogonal to the hypersurface. The vector 
n is collinear with the timelike normal evolution vector \/t as follows: 



whereas the 1-form, n, which is dual to n, is collinear with the gradient 1-form dt, 
as follows: 



n = v t, 



(2.1) 



n = -Ndt. 



(2.2) 



N is thus the lapse function that ensures the normalization relation: 



n ■ n = (n, n) 



= -1. 



(2.3) 
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On each hypersurface, we introduce a spatial coordinate system ( 
which constitutes a well-behaved 4- dimensional coordinate system on M, (x a ) = 
(t, x 1 , x 2 , x 3 ) when varied smoothly between neighboring hypersurfaces. The natural 
basis for this coordinate system is denoted by (d a ) = {d t , <%) where: 

9t '= Ft 

di := ^,16 1,2,3. (2.4) 

d t is the time vector tangent to the lines of constant spatial coordinates, whereas 
di,i e 1,2,3 are vectors tangent to the hypersurface. We define as "coordinate 
observers" a class of observers whose 4-velocity is collinear with the time vector. 

The difference between the time vector dt and the timelike normal evolution vector 
\/t is the shift vector f3, and the relation is as follows: 

dt = ^t + (3 = Nn + /3. (2.5) 

The shift vector thus enables the freedom to choose how the spatial coordinate system 
changes from hypersurface to hypersurface. Fig. 2.1 shows the construction of an 
adapted coordinate system. 

The 4-metric g can be written in terms of components g a p with respect to the 
coordinates (x a ) as follows: 

g a p = g(d a ,dp). (2.6) 



Using Eq. (2.5), the 00-component of the 4-metric is thus given by: 



g 00 = g(0 t , d t ) = dfd t = -N 2 (n-n)+(3-(3 = -N 2 + (3 t (J\ (2.7) 
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whilst the Oi-components are given by: 

9oi = g(d t , d t ) = (Nn + /3)-di = ft, (2.8) 

as the vector di is tangent to the hypersurface. The 3-metric 7 is thus induced on 
each hypersurface with respect to the adapted coordinate system via this relation: 

Qij = g(0i, 8j) = j(d h dj) = 7^. (2.9) 

The orthogonal projector that maps the 4-dimensional metric g onto the 3-dimensional 
hypersurfaces is given in terms of components with respect to the coordinates (x a ) 
by: 

tf = 6% + n a nf,. (2.10) 
The operator acts on the normal timelike unit vector, as follows: 

7 «n" = 5^n p + n a n p nP = n a - n a = 0, (2.11) 

and on any vector tangent to the hypersurfaces, as follows: 

tfd a = 6p a + n a n fi d a = df ) . (2.12) 

2.2 The covariant and Lie derivatives 

The 4-dimensional spacetime M with the metric g possesses an associated connec- 
tion, V) denoted as the affine connection or covariant derivative, which enables the 
comparison between vectors evaluated at 2 different points along the congruence of 
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curves generated by a vector field. We define the vector evaluated at a point P, 
with coordinates x a along the congruence, as Vp and the vector evaluated at a point 
Q, with coordinates x a + 5x a , as Vq. This vector at point Q is given by a Taylor's 
expansion as follows: 

v a (x + 8x) = v a (x) + 8x p d p v a . (2.13) 

We further introduce another vector at point Q which is 'parallel' to the vector at 
point P, and denote it as Vp + 5v a . 5v a is collinear with Vp and 5x a and thus can be 
written as: 

5v a (x) = -T a ^{x)5x ri . (2.14) 

The connection or covariant derivative then evaluates the difference between vf> and 
Vq as follows: 

Va^ = -L K - v%] = v ? a + r (2.15) 

where are thus denoted as the connection coefficients, also known as Christoffel 
symbols. The covariant derivative can be generalized to apply to the differentiation 
of a tensor T of any rank in the 4-dimensional spacetime M, along a congruence 
generated by any vector field u. The 3-dimensional covariant derivative, DT, is thus 
obtained from the projection of the 4-dimensional covariant derivative, on to the 
3-dimensional hyper surf aces, as follows: 

DT = 7VT. (2.16) 

The connection coefficients components can be given in terms of the 4-metric 
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components as follows: 

T lf3 = \9^[9an,P + 9^,a ~ 9«pJ- (2.17) 

This can be written analogously for the 3-dimensional connection coefficients in terms 
of the 3- metric components, 7^. 

The 4-dimensional Lie derivative for M measures the distortion of the 4-dimensional 
coordinate system. As it is based solely on coordinate bases, the 4-dimensional Lie 
derivative is considered a more fundamental construct compared to the affine connec- 
tion. The Lie derivative is given by the difference between the vector evaluated at 
point Q, Vq, and the vector evaluated at point P that is 'dragged along' to point Q, 
given by: 

v%(x + 8x) = v' a (x + 5x) 

= v' a (x + 5su a ) 
dx' a a 

= V 

dx? 

= {51 + 6sd p u a )v p {x) 

= v a (x) +5s(d u a vP(x), (2.18) 
where s is the affine parameter along the congruence of curves generated by the vector 
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field u. Therefore, the derivative is given by: 

C u v a = lim ^-[v a (x + Sx) - v' a (x + 5x)} 

Ss^O OS 

= lim -j-[v a (x) + dsu^daV 01 - v a (x) - 5s(9„-u Q -u M (a;)] 

Ss^O OS 

= u^v -d^v". (2.19) 

The 3- dimensional Lie derivative acts on vectors tangent to the hypersurfaces in the 
same way. 



2.3 The intrinsic and extrinsic curvatures 

The intrinsic curvature of the 4-dimensional spacetime M is described by the non- 
vanishing of the commutator of any vector, v in M, (y a V/3 ~~ V/3 Va)^ as follows: 

(Va V/3 - V/3 V>" = V[a V/3] ^ = fl?a^ 7 . ( 2 - 20 ) 

Again, this can be generalized to a tensor of any rank in M. The contraction of the 
intrinsic curvature, also known as the Riemann tensor R^ aj3 , gives the Ricci tensor 
FLy/3. The contraction of the Ricci tensor in turn gives the Ricci scalar R = g" //3 R 7 /s. 
This 4-dimensional Ricci scalar is independent of the ambient coordinate system of M. 
The expressions for the 3-dimensional Riemann tensor, Ricci tensor and Ricci scalar 
take on analogous forms. We shall denote them as ^R\j k , ^Rij and ^R respectively. 
Similarly, the 3-dimensional Ricci scalar is independent of how the hypersurfaces are 
embedded in M. 

17 



Chapter 2 3+1 formalism in numerical relativity 

Conversely, the extrinsic curvature describes how the hypersurfaces are embedded 
in M, or more specifically, it measures the curvature of the hypersurfaces in the 
embedding. It is given by the orthogonal projection of the covariant derivative of the 
timelike normal unit vector n along any vector u tangent to the hypersurfaces, as 
follows: 

K = - 7Vu n. (2.21) 
By invoking the vanishing of the covariant derivative of the 4-metric as a direct result 



of Eq. (2.17), the covariant expression of the projection operator in Eq. (2.10), and 



the n a \j u n a = identity, the extrinsic curvature can also be expressed as the Lie 
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derivative of the 3-metric along n as follows: 

= (( V „JV) Vat + N V „ Va t)^ 

= D p N\7 a t + N\/ a (-N~ l n^ 

= -N-WDpN - y a {N-\Nn^) - \J a n^ 

= -n a D p In N - (v a nji{ - ( XJ (an^np 

= -n a D/3 In N - \J a np 

= - g k"" ™/3 + ™/3« M Vm "a + V« n /J + S7pn a ] 

= ~^[n a n^ v M + npn^ Vm "a - n a n^ V/3 % ~ n /3 n ^ V« »fy + Va(ff?X) + 

= - 2 K (w a n/j) - n a n M V/3 % - n^n" V« »fy + Qfr V« ™ M + 9an V/3 ™1 
= (fiW* + ^ n /3) + (gpn + npn^) Va n' 1 + (g afM + n a n^) V/3 n^] 

= ~ ^ K 7a/3 + 7/3 M Va ^ + 7a M V/3 ™1 

= —Cny. (2.22) 
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2.4 The Gauss-Codazzi relations 

The various contractions of the 4-dimensional Riemann tensor ( 4 )R along the 3- 
dimensional hypersurfaces and the timelike unit normal vector n are essential for 
the 3+1 formulation of the Einstein field equations. Solutions to the field equations 
can be obtained when the equations are cast in an initial value problem involving 
constraint equations on the initial hypersurface and evolution equations on subse- 
quent hypersurfaces. Due to the fact that these contractions of ( 4 )R do not involve 
any timelike derivatives of the metric tensor g, they will be used to construct the 
constraint equations. In this section, we present how these contractions are obtained. 
Covariant derivatives of vectors and tensors will hereby be denoted by the subscripted 
semicolon, and the vectors and tensors along the hypersurface will use indices denoted 
with the alphabet instead of Greek letters. 

We begin by considering the full projection of onto a hypersurface. Using 



Eq. (2.10) and Eq. (2.21), we first calculate the commutator of the double covariant 
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derivative of the projection tensor 7 as follows: 

7^7? 7 C 7 ~ ll^llll = -T d bc T e ad i: + T d cb T e adl : + Y d bc K ad n a - Y d b K ad n a + 

r^(r^ e a ) - r d ac (r dbl :) - r d ab K dc n a + r d ac K db n a - 

K ab , c n a + K aCyb n a - K ab n^ c + K ac n^ 

/■pe pe 1 -pcf pe pd pe \~/Ot _ 

V afc,c ac,6 a& dc ac dfe/ le 

(K ab , c - T d ac K db )n a + (K ac<b - Y d ac K dc )n a - 
K ab n^ c + K ac n^ b 
= -K bc < - (K ablc - K ac \ b )n» - 

K ab n^ c +K ac n^. (2.23) 

Denoting analogously to ^-R^,, the 3-dimensional intrinsic curvature ^R]j k , in 
terms of vectors v l tangent to the hypersurface as follows: 

-•"V^iri ( 2 - 24 ) 
the previous equation can be written as: 

(4) K^:ibl2 = (3) K bc < + {K ab , c - K ac , b )n» + K ab n^ c - K ac n^ b . (2.25) 

A further projection of this equation along the hypersurface yields the Gauss relation 
as follows: 

{i) K^a^H = (3) RL + (K d K ac - K d K ab ). (2.26) 
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The Gauss relation can be contracted in the /i and /3 indices using again Eq. (2.10) 
to yield the following: 

7a a 7c ^i4 7 + 7cX7M 4) ^/3 7 = (3) Rac + KK ac - K\K ah . (2.27) 

Further contraction in the a and c indices will yield the contracted Gauss relation as 
follows: 

(4) i? + 2 (4) i? Ma n"n a = (3) R + K 2 - K ab K a \ (2.28) 
which will be used in the 3+1 formulation of the Einstein field equations. 



We now consider the projection of Eq. (2.26) along the timelike normal unit vector 
n which results in the Codazzi relation: 

^KfnnrfTb 7c 7 = K ab-,c - K ac , b . (2.29) 

The Codazzi relation can be similarly contracted in the a and b indices to yield the 
contracted Codazzi relation: 

= K. c - K* (2.30) 



2.5 The Ricci relation 

The projection of ^ 4 ^R twice along the hypersurface and twice along the timelike unit 
normal vector n yields the Ricci relation which will be used to construct the evolution 
equations in the 3+1 formulation of the Einstein field equations. We obtain this by 
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considering the projection of the commutator of the double covariant derivative of the 



timelike unit normal vector itself, and using the expression for K a p in Eq. (2.22) and 



the component form of Eq. (2.21), as well as the projection onto the hypersurface of 
the Lie derivative of K a p with respect to the timelike normal evolution vector m, as 
follows: 

= 7«Xt£[- W + DnnNn a ), v + (K£ + D»lnNn v ), a ] 
= ^n a ^[-K-,u - n a , v DHnN - n a {mnN). u - - 

n v . a DHnN - n u (DnnN), a ] 
= l a ^[-n a K^ u + {D»lnN). v + n°K^ a + n°n U]a D»lnN] 
= lc&ifa K >% + [DnnN), v + n°K». a + D v lnNDHnN] 
= -K aa Kl + DpDJnN + ^rfK^ + DJnNDplnN] 
= -K aa Kl + ^DpD a N + i^tfK^ 
= -K aa K% + -j^DpDaN + ^{C m K aP + 2K aii K%) 
= -j^DpDaN + -j^C m K al3 + K aix Kp. (2.31) 



2.6 Projections of the stress-energy tensor 

The various projections of the stress-energy tensor onto the 3-dimensional hypersur- 
face and along the timelike unit normal vector n are needed to construct the matter 
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sources for the 3+1 formulation of the Einstein field equations. 

We first present the full projection of this tensor along n. We recall from Section 
2.1 that the 4- velocity of the Eulerian observers is definitionally the timelike normal 
unit vector n. Therefore, the projection of the stress-energy tensor along n is the 
matter energy density, which is a scalar measured by these Eulerian observers. We 
denote this matter energy density as E. 

The mixed projection of the stress-energy tensor is called the matter momentum 
density, p, which is a linear form tangent to the hypersurface. Its component form is 
given by: 

Pa = -T^n u . (2.32) 

The full projection of the stress-energy tensor along the hypersurface is called the 
matter stress tensor, S, which is a bilinear form tangent to the hypersurface. Its 
componenet form is given by: 

Sap = T^lalp- (2-33) 



Using Eq. (2.10) in this component form of the matter stress tensor and taking 



the trace, we obtain the following relation: 

T = S-E. (2.34) 
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2.7 3+1 decomposition of the Einstein field equa- 
tions 

Armed with the various projections of the intrinsic curvature tensor ^R and the 
stress-energy tensor S, we are now ready to present the full 3+1 decomposition of the 
Einstein field equations. In the decomposition, we will utilize two forms of the field 
equations, namely: 

1^ 

{4) R - - Rg = 8ttT, (2.35) 
2 



and its equivalent: 



^R = 8n(T-^T g ), (2.36) 



where T := g^T^ is the trace of the stress-energy tensor T. 

To construct the Hamiltonian constraint components of the Einstein field equa- 
tions, we apply the twice-contracted Gauss relation obtained in Section 2.4, into the 



twice-contracted Eq. (2.35) along the timelike normal unit vector n, as follows: 



1(4) 

R^nW - - Rg^n v = ZitT^n" 

2R^n»n u + (4) R = 16nE 

{3) R + K 2 - K ab K ab = IQttE. (2.37) 

The momentum constraint components of the field equations however will make 
use of the contracted Codazzi relation as obtained in Section 2.4, in the mixed pro- 
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jection of Eq. (2.35) once along n and once along the hyper surf ace, as follows: 



1(4) 



8np c . 



(2.38) 



To construct the evolution components of the field equations, we first combine the 
Ricci relation (Section 2.5) with the once-contracted Gauss relation (Section 2.4): 



1^ {4) R^ = (3) R af s + KK afS - -^DpD a N - ^C m K a/3 - 2K afM K%. 



N 



(2.39) 



We then substitute this into the equivalent form of the field equations Eq. (2.36) in 



this construction: 



^R a p + KK a/} - ^D p D a N 



N 



C m K a p-2K a „K$ = 8 7r [S a0 --(S-E) laP ] 

£ m K a p = —DpD a N + N{^R a/3 + KK a/3 



2K afl K% + 4vr[(S - E) la p - 2S Q/3 ]}.(2.40) 
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3.1 3+1 decomposition of general relativistic hy- 
drodynamics 

The conservation of the stress-energy tensor T ensures that, as long as they are 
satisfied on the initial hypersurface, the Einstein field equations are satisfied for all 
time. The conservation is given as: 

Tg = 0. (3.1) 

In neutron star simulations, the neutron star matter is characterized by a perfect 
fluid, where there is no presence of heat conduction and other stresses besides pressure. 
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The stress-energy tensor for a perfect fluid is purely defined by the matter energy 
density and the pressure. In component form, the stress-energy tensor is given by: 



T/xi/ = phu^u u + pg^. 



(3.2) 



To incorporate these equations into the 3+1 formalism, we decompose them into 
the 3+1 form by first writing them in a 1st order flux conservative form: 



d t U + 8iF l = s, 



(3.3) 



where U, F l and S are the evolved state vector, the flux vector and source vector 
respectively. The evolved state vector U is written in terms of the primitive variables, 
ie. the matter density p, fluid velocity vector v\ and specific internal energy density 
e as follows: 



(3.4) 





D 




^fWp 


u = 






yfyphW 2 Vj 




T 




y/l(phW 2 -P- Wp) 



where the fluid velocity vector is a 3-velocity related to the 4-velocity u l as: 



W 



m = jj{i,Nv*-f?}, 



(3.5) 



with W representing the Lorentz factor W — 1/ a/1 — %jV i v : >, and h is the specific 
enthalpy which can be written in terms of the specific internal energy density as 
follows: 

h = 1 + e + P/p. (3.6) 
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The flux vector F % is defined as: 



N(v* - fi/N)D 



(3.7) 



7V((V - ^/N)t + v/tVP) 



whilst the source vector S is defined as: 







S 



(3.8) 



3.2 Conformal decomposition 

Conformal decomposition was introduced by Lichnerowicz in 1944 [33] to facilitate a 
more efficient resolution of the constraint equations in obtaining valid initial data for 
the initial value problem. In the decomposition, the 3-metric 7 is written in terms of 
a conformal factor which is a positive scalar field, and a conformal 3-metric 7 as 
follows: 



where det 7 = det / = 1 when Cartesian coordinates are used. The extrinsic curvature 
K of the 3-dimensional hyper surf ace is decomposed into its trace K := K\ = 
and a traceless form A as follows: 



7 = ^ 4 7, 



(3.9) 




(3.10) 
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where 7^A« = 0. We recall that the expression of the extrinsic curvature can be 



written in terms of the Lie-derivative of the 3-metric, as shown in Eq. (2.22). 



We substitute the new form of K as written in Eq. (3.10) and 7 into this to obtain 



the following: 



^m(^ 7ij) 



(3.11) 



Taking the trace of this equation with respect to 7 and applying the general law of 
variation of the determinant of an invertible matrix twice, we obtain the evolution 
equation for the conformal factor \I/ as follows: 



^-£ fl )ln*=±(D i p i -NK). 



(3.12) 



Using the general law of variation of any invertible matrix, Eq. (3.11) also yields 
the evolution equation for the conformal metric 7 as follows: 



,d_ 
dt 



-2iV^- 4 A 



2 ~ 



(3.13) 



In order to obtain the conformally-decomposed evolution equation for the extrin- 
sic curvature tensor K, we recall the evolution component of the Einstein field equa- 



tions in the 3+1 formalism as obtained in Chapter 2, ie. Eq. (2.40), and substitute 



Eq. (3.10) into it. Also using Eq. (2.22), we obtain: 



2 

C m Kij = C m Aij + -C m K^ij - -NKKij. 



(3.14) 



30 



Chapter 3 



Application of the 3+1 formalism in neutron star simulations 



From the identity: 



fi£ m K ij + 2NK ij K ij , 



(3.15) 



we thus obtain: 



C m K = —D i D i N + N[R + K 2 + Att(S — 3E)] 



-DiD l N + N[4tt(E + S) + K 



(3.16) 



We now substitute this and Eq. (2.40) back into Eq. (3.14), employing again Eq. (3.10) 
to obtain: 



r a .. 



1 



1 



-DiDjN + N[Ri3 + -KA {j - 2A lk A) - 8vr(^ - -S 7ij )) + 



\{D k D k N -NR) liy 



(3.17) 



Eq.s (3.16) and (3.17) represent the trace part and the traceless part of the evolution 



equation for K. We now conformally decompose the trace part, ie. Eq. (3.16) by 



subtituting Eq. (3.10) and A %1 = ^ A %3 into it to obtain the following: 



_ _ _ _ K 2 

C m K= -m A (DiD l N + 2DMmD l N) + N[Aix(E + S) + A i5 A lj + — ]. (3.18) 

o 



Similarly, employing Eq. (3.12), the conformal connection C t * = — \l k \-^r + 



dxi dx 



x) and the resulting conformal Ricci tensor and conformal Ricci scalar, we 
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conformally decompose the traceless part into: 



r 4 



—D^Aij + N[KAij - 2f l A ik A jt - 8tt(^-% - ^S%)} + 

ty~ 4 {—DibjN + 2Di In VDjN + 2D d In VDiN + 

\(b k D k N - AD k In *D k N)% + 

N[Rij - -Rlij - 2DiDj In * + 4A In VDj In * + 
3 

§ In * - 2.Dfc In ^D fe In } . (3.19) 



We now present the conformal decomposition of the constraint part of the Ein- 



stein field equations in the 3+1 formalism as obtained in Chapter 2, ie. Eq.s (2.37) 



and (2.38). We substitute Eq. (3.10) and the conformal Ricci scalar, R = \& 4 R — 
8ty~ 5 DiD l ty, into the 3+1 Hamiltonian constraint equation, ie. Eq. (2.37), to yield: 



1 s 



1 7 



1 



DiD l V - -RV + (-Ah A* 3 K 2 + 2nE)^ b = 

8 V 8 3 12 ; 



(3.20) 



Similarly, the conformal decomposition of the momentum constraint equation, ie. 



Eq. (2.38), yields: 



DjAWj In * - -&K = 8tt* V- 



(3.21; 



These six equations, ie. Eq.s (3.12), (3.13), (3.16), (3.19), (3.20) and (3.21), con- 



stitute the conformal 3+1 Einstein field equations. This set of equations is solved 
for the conformal 3-metric 7^, the conformal traceless part of the extrinsic curvature 



Aij, the conformal factor \l/ and the trace of the extrinsic curvature K. We recover 
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the physical 3-metric 7 and the physical extrinsic curvature K via the following: 

lij = ^ lij 

K l3 = ^\A l] + l -K^). (3.22) 

However, in our neutron star simulations, we employ a modified form of the 
above conformal 3+1 Einstein field equations. This modified formulation is called 
the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation. In this formulation, 
a vector is introduced by Shibata and Nakamura as well as Baumgarte and Shapiro 
to restore the Laplacian nature of the conformal Ricci tensor that is written in terms 
of the conformal metric. 

In order to do this, the expression of the conformal Ricci tensor in terms of the 
conformal metric, 7, is considered as follows: 

f) _ f) _ _ _ _ _ 

p _ ~pk _ -pk _|_ -pfc -pi -pk-pl (o oo\ 

j ~ dx k ij dxi ik v kl il k i' ^ ' 

We introduce a new tensor field A as follows: 



A* := f* - f*, (3.24) 



where T k j are the connection coefficients for the flat metric with respect to the coor- 
dinates (x % ). The components of this tensor field can also be written as: 

A* = \l kl m l3 + Vfi u - 2>,7y), (3.25) 

where T>i represents the covariant derivative associated with the flat metric. 
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Using the expression (3.24) together with the fact that the Ricci tensor vanishes 



for a flat metric and the identity A k k = 0, we obtain: 



^ ~ L A * fc + A ^ AL + ***** + ^ AL " ^ 

pi A k pfc A I 

1 kj^il 1 il^kj 




A iA l kj 



Afc I pi Afc pi Afc P^A' 

Q x k V + 1 kl^ij 1 kj^il 1 il^kj 



(3.26) 



Using expression (3.25) on the previous result yields: 



Rjj 



U) k \f\Va i3 + Vfju - Vjriij)] - AjA^ 

v^ay - f l v k va l3 } - AX, 



zM 



zM 



zM 



vwfvfm - f l v k v^) - a* Al, 
l 



il kj 
~,kl i ~ ~,kl\ 



^ (TVkDrtij + likVpiT + i Jk v t Va kl ) + o fi ( 7 , x> 7 ), (3.27) 

where (^(7, D7) := -i(D fc7 ^7 fci + V^uVrf 1 + V^Vflij) ~ A£Ajy. 

The Ricci tensor as rendered in the previous expression can be viewed as a Laplace 
operator acting on the conformal metric 7 yielding second-derivatives on the right 
hand side. However, the second and third terms on the previous expression spoils 
the elliptic character of the Ricci tensor operator. It is here that Baumgarte and 
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Shapiro introduce the vector P = —Dff^i which turns the previous expression into 
the following: 

Rij = l -(-i kl V k Va l3 + 7^f fc + 7 jfc Z\f k ) + Oijtf, TVy). (3.28) 

Taking the trace of this expression of the conformal Ricci tensor, as well as recalling 
the identity j^'Dijij = 2Af fe = 0, the conformal Ricci scalar is thus: 

R = V k f k + 0(j,Vj), (3.29) 

where 0(7, T>j) := \^V k ^T>i^ij + 7^0^(7, T>j), which is a term that does not 
contain any second derivatives of 7 and is quadratic in the first derivatives. 

Earlier, Shibata and Nakamura introduced the covector Fi = T>^ij instead of the 
vector P, which is related to the latter as follows: 

Fi = %-P - (f k - P k )V k ^ r (3.30) 

However, the vector P has an edge over the covector Fi because it covers all the 
second derivatives of the conformal metric that do not contribute to the Laplacian 
operator. 

3.3 Gauge choices 

Looking at the set of conformally decomposed Einstein field equations obtained in 
the previous section, we note that there are no derivatives of either the lapse function 
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N or shift vector (3. This indicates that in the solution of the field equations, the 
lapse function and the shift vector can be freely chosen while yielding the same phys- 
ical solution. We recall from Chapter 2 that the lapse function determines how the 
spacetime is sliced and the shift vector determines the choice of coordinates on these 
spacetime slices. The choice of the lapse and the shift thus changes the form of the 
Einstein field equations to be solved, making it either more hyperbolic or more elliptic 
in nature. As the success of numerical simulations in modeling neutron star systems 
depends crucially on the well-behavedness or non-singularity of the coordinate func- 
tions, the freedom of the lapse and shift choice gives us the privilege to adjust the 
hyperbolicity of the field equations based on the nature of the physical system to 
be studied. We will discuss in this section the common choices made particularly in 
neutron star simulations. 

The simplest lapse choice is called the geodesic slicing, which sets N — 1. Since 
the acceleration co-vector can be given by: 

a a = D a ]nN, (3.31) 

setting N — 1 renders zero acceleration for the Eulerian observers, ie. they travel 
along geodesies of the spacetime, hence the name of this lapse choice. This choice 
permits only limited evolution of the spacetime due to its focusing property that 
results in coordinate singularities. 

A more popular choice is the maximal slicing which sets the extrinsic curvature 
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scalar K = 0. This lapse choice maximizes the volume of the hypersurface. The 
volume enclosed within a closed 2-dimensional surface lying on a hypersurface is 
given by: 

V = I ^d 3 x, (3.32) 
Jv 

with 7 is the determinant of the metric 7 with respect to the coordinates (x*) on the 
hypersurface. The change of hypervolume can thus be written: 



dV f d^f 



dt J Vt dt 



d 6 x, (3.33) 



where Vt is the domain at time t. Using the identity 7^~g^ = —NK + Di/3 1 that is 
derived from the evolution equation of the 3-metric, the variation of the hypervolume 
in the previous equation can be written as: 



f [-NK + D^]^d 3 x 

JVt 

NK^d 3 x + j> (3 l Si^qd 2 y 



dV 
~dt 



NK^d 3 x, (3.34) 

v t 

where s is the unit normal to S lying in the hypersurface, q is the induced metric 
on S with q = q a b and (y a ) are the coordinates on S. Therefore, setting K = on 
the hypersurface renders the hypervolume extremal with respect to the variations in 
the domain bounded by S. With a metric of a Lorentzian signature, this extremum 
is a maximum, hence the name of this lapse choice. Combining the maximal slicing 



condition with the evolution equation for the extrinsic curvature scalar K, Eq. (3.16) 
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as derived in the Section 3.2, we obtain: 



DiD l N = N[An(E + S) + h'.jh'' 1 



(3.35) 



an elliptic equation which imposes the condition on all subsequent hypersurfaces in 
the spacetime. Maximal slicing avoids the formation of coordinate singularities in the 
evolution of a physical system. As we recall from Chapter 2, the extrinsic curvature 
scalar is defined as the divergence of the unit normal vector n, alternatively called the 
4- velocity field of the Eulerian observers. Fixing K = thus results in y-n = 0, which 
prevents the Eulerian observers from converging towards any coordinate singularity 
that forms due to the focusing effect of gravity. This happens when the lapse function 
and thus the proper time between two adjacent hypersurfaces tends to zero as the 
coordinate time tends to infinity. 

The next important category of lapse choice is called the 1 + log slicing, which 
is a generalization of the harmonic slicing introduced by Bona, Masso, Seidel and 
Stela [Hj. Harmonic slicing was introduced by Choquet-Bruhat and Ruggeri [2T] 
in an attempt to write the 3+1 Einstein field equations in a hyperbolic form, and 
sets Dgt = 0, a harmonic condition for the time coordinate. Using the relation 
y/— g = Ny/j, the expression for the 4-metric g al3 delineated in Chapter 2 and again 
the identity -^=-^r = —NK + Di(3 l , this d'Alembertian becomes: 



which is an evolution equation for the lapse function. The 1 + log slicing thus gener- 



di 



Cp)N = -KN\ 



(3.36) 
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alizes this equation as follows: 

(^--Cp)N = -KN 2 f(N), (3.37) 

where / is an arbitrary function with f(N) = 1 corresponding to the harmonic slicing 
and f(N) = corresponding to the geodesic slicing. Using the identity — 
—NK + DiP 1 and employing normal coordinates where f3 — 0, this equation becomes: 

dN d , , 

«■ = «'■"'• (X38) 

which yields N = 1 + In 7 as one of its solutions, hence the name of 1 + log slicing. 
The 1 + log slicing produces foliations very similar the to maximal slicing and thus 
have strong singularity avoidance, its biggest advantage, making it another popular 
choice for neutron star simulations. 

We now move on to the gauge choices commonly made in the shift vector for neu- 
tron star simulations. Again, the simplest of these choices are the normal coordinates, 
which set (3 = 0, where the 4- velocity field n or the unit normal vector field lines are 
parallel to the constant spatial coordinate field lines, hence the name normal coordi- 
nates. An advantage of this choice is its incapability in introducing any pathologies 
of its own. However, in rotating star spacetimes, the employment of this shift choice 
can result in coordinate shears due to the fact that the unit normal vector field lines 
are not parallel to the stationary Killing vector field £ lines. 

The minimal distortion shift was introduced by Smarr and York in 1978 [75] in an 
attempt to minimize the time derivative of the conformal 3-metric, 7 := Cq^. The 
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components of this time derivative are: 

1% = % (3-39) 
which are related to the distortion tensor components Qij as follows: 

Qa = * 4 ^, (3-40) 



which has 5 degrees of freedom, ie. 6 freely chosen components minus 1 constraint 
requiring det 7^ = 0. The distortion tensor Q is a measure of the change of the shape 
of spatial domain V within a fixed coordinate boundary from one hypersurface to 
the next. The minimal distortion shift hence seeks to minimize the change in the 
shape of the spatial domain V. This can be done by choosing coordinates (x l ) that 
set the distortion tensor Q identically equal to zero. Taking into account that the 3 
degrees of freedom for the coordinate choice are insufficient to constrain the 5 degrees 
of freedom for the distortion tensor Q, we decompose the distortion tensor Q into a 
longitudinal part and a transverse-traceless part as follows: 

Q tl = (CX) lJ+ Q^, (3.41) 

where CX is the conformal Killing operator associated with the physical metric 7 
acting on some vector field X. Taking the divergence of this relation, we obtain: 

I )>(.),, D'iCX);,. (3.42) 

which possesses 3 degrees of freedom that can be constrained by the coordinate choice. 
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Hence, the minimum distortion shift conditon becomes: 



D j Qij = 0. (3.43) 



Using the evolution equation for the 3-metric 7, namely, Cm^ij = —2NKij, the 
identity ^¥ = ~ NK + and the 

expression for the traceless part of the 



extrinsic curvature tensor A, Eq. (3.10), this condition yields: 



- 2NDjA ij - 2A ij D 3 N + + \ - 2D k fi k lid ) = 0. (3.44) 

Furthermore, we can employ the 3+1 momentum constraint equation, Eq. (2.38), to 
obtain the following elliptic equation of shift evolution resulting from the minimal 
distortion condition: 

DjDlfr + X -D L DSi + m/3 j = 16rrNp i + -ND { K + 2A ij DjN. (3.45) 
3 3 

A more computationally-feasible shift choice based on the minimal distortion shift 
is the T-freezing shift, introduced by Alcubierre and Briigmann [3]. This shift choice 
sets Dj7 iJ = -^(Vj'j^) = 0. The covariant derivative Dj can be written in terms of 
the connection coefficients for the flat metric T l - k with respect to the coordinates {x l ) 
as follows: 

Vfp = f k (T) k - f) k ) = -r. (3.46) 

Hence, the T-freezing shift condition sets = 0. To apply this condition on the 
shift evolution, we write the Lie derivative of the conformal 3-metric 7^ in terms of 
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the covariant derivative as follows: 

ft = 2N& + f3 k V k f j - f j V k /3 k - f k V k ft + 2 -V k (3 k f\ (3.47) 

the divergence of which, via the expression T>jfi = — P, the conformal 3+1 mo- 
mentum constraint equation, and the relationship between the conformal and flat 
covariant derivatives, becomes: 

dt 1 ~ ~ 2~ 

— = -2NVjA tJ - 2A ij VjN + (3 k V k T l - T k V k (3 l + -T l V k f3 k + 

f k vp k p + l -f 3 vp k $ k , 
f k vp k p + \fwp k (3 k + \rv k p k - Y k V k fi + (3 k V k r 

3 3 
= 2iV[87r*V - ^' fc (r}fc - T) k ) - 6A j Vj In * + \f 3 VjK\ + 

2A ij VjN. (3.48) 

The T-freezing shift condition thus yields an elliptic equation for the shift. Alcubierre 
and Briigmann [3] turned it into a parabolic equation using the relation ^j- = k^- 
with k being a positive constant. A modification of this is used in the neutron star 
simulations presented in this thesis, namely: 

JL = df - C 2 (3\ (3.49) 

where we set C\ — C2 — 0. 
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3.4 The initial value problem 

In Chapter 2, we see that the Einstein field equations are separated into the constraint 
part and the evolution part, such that their resolution amounts to solving an initial 
value problem using the constraint part to obtain initial data that will be propagated 
forward in time using the evolution part of the field equations. As the initial data is 
constrained, it is a non-trivial astrophysical problem common particularly in neutron 
star simulations, to ascertain that the solution to the constraint part of the field 
equations yields the physical system to be studied. 

In this section, we discuss the conformal transverse traceless method proposed 
by York [83j,[84j,[85j, a method that has been employed in our general relativistic 
hydrodynamics simulations of neutron star and neutron star-like systems. We recall 
that in Section 3.2, we have obtained the conformally decomposed constraint part of 
the Einstein field equations. In 1973 and 1979, York solved this system of equations 
by further decomposing the trace part of the conformal extrinsic curvature A % i into 
a longitudinal part and a transverse part, as follows: 

A ij = (CX) lJ + A% T , (3.50) 

where A^ T is both traceless and transverse with respect to the conformal metric, ie. 
jijAj, T = and DjA^ T = respectively, and (CX)^ is the conformal Killing operator 
associated with the conformal metric 7 that acts on the vector field X as follows: 

{CX) ij := &X j + D j X i - -D k X k f j , (3.51) 

3 
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and is also traceless, ie. ^^(CX)^ = 0. 

We incorporate this longitudinal/transverse decomposition and the identity A^X* 
DjA 1 ^ into the conformally decomposed constraint part of the Einstein field equations, 



Eq.s (3.20) and (3.21 ) and obtain the following: 



D i b i ^-\m+\[{CX) ij +AjT][(CX) i i+ = 0, (3.52) 



A c X l - -y 6 &K = 8np l 
3 



where (CX)^ := 'y ik 'y jl (CX) kl and A\( = ^ ik 1jiA^ T . 



(3.53) 



Eq.s (3.52) and (3.53) above show that in solving this system of equations, the 



conformal metric 7, the symmetric traceless and transverse tensor Af T , the extrinsic 
curvature scalar K, and the conformal matter variables (E,p l ), can be freely chosen, 
whilst the conformal factor \1/ and the vector X will be determined as results of the 
initial value solve. We then construct the physical metric as 7^ = the physical 

extrinsic curvature tensor as = \I/ -10 ((£X) Jjf + A^ T ) + ^~ 4 K^ , the physical 
matter energy density as E = ty~ 8 E and the physical matter momentum density 
vector as p % = ^/~ 10 p\ which form a set of initial data (7,K,i?,p) that satisfies the 



constraint part of the Einstein field equations, Eq.s (2.37) and (2.38). 
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3.5 GRAstro-2D as an axisymmetric general rela- 
tivistic hydrodynamics code 

The GRAstro-2D code is based on the Cactus Computational Toolkit [UJ and the 
GRAstro-3D code [3"T] . |3"2"] . Thus, similar to the GRAstro-3D code, it solves the full 
3+1 Einstein field equations as presented in Chapter 2, using the BSSN scheme as 
presented earlier in Section 3.2. The evolution of the code is similarly unconstrained, 
ie. the constraint equations are only solved at the initial time and not throughout the 
evolution. As such, violation of the constraint equations are monitored throughout 
the evolution to determine convergence of the code. The code also employs the same 
initial value problem solver and finite-differencing scheme as used in GRAstro-3D. In 
this section, we present the basic concepts behind the modifications made by [UJ on 
the GRAstro-3D code for the purpose of adapting the former to perform numerical 
simulations of axisymmetric systems. In both GRAstro-3D and GRAstro-2D, we use 
geometric units where we set G = c = 1 (refer to Appendix C). 

The basic concepts are based on the Cartoon technique introduced by Alcubierre, 
Brandt, Briigmann, Holz, Seidel, Takahashi and Thornbug in 2005 [4J. This technique 
is able to avoid the problems normally encountered in other axisymmetric techniques, 
eg. techniques that employ a cylindrical (p, z, 9) coordinate grids, where physically 
non-singular variables may become indeterminate 0/0 forms along the z-axis. Al- 
though these indeterminate forms can be regularized by L'Hospital's rule, their eval- 
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Figure 3.1: Figure taken from [3]. 



uation becomes inaccurate in finite differencing schemes where the grid spacing is 
finite when its limit is required to go to zero to be consistent with its analytic coun- 
terpart. In addition to this problem, some of the variables in the 3+1 Einstein field 
equations are obtained in finite differencing schemes using a summation of terms 
that may include these indeterminate forms. Regularization of such variables require 
detailed analysis of the entire system of the 3+1 Einstein field equations in the ax- 
isymmetric coordinate system near the axis of symmetry, a difficult although not 
impossible undertaking. 

However, Cartesian coordinate grids do not introduce any such pathology, as the 
coordinate system is completely regular even at the grid origin. The Cartoon tech- 
nique employs a 3-dimensional Cartesian (x, y, z) coordinate grid that is restricted to 
the y = plane by only one finite-diffence-length. It uses continuous rotational sym- 
metry to provide the boundary conditions for this thin 3-dimensional slab. Fig. 3.1 
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shows the construction of this slab, where the z axis is taken as the axis of symmetry. 



The solid dots represent grid points where variables take on the values calculated 
on the original 3-dimensional coordinate grid. The white dots represent grid points 
where the variables are calculated via the Cartoon boundary conditions. The white 
squares represent grid points where the variables are calculated by imposing a physical 
boundary. The Cartoon boundary condition is described by a rotational coordinate 
transformation R. If we consider arbitrary tensor fields T on the 3-dimensional co- 
ordinate grid, a rotation of such a tensor field by an angle — 6q about the axis of 
symmetry is equivalent to the rotation of the coordinate system by an angle 8 . The 



coordinate transformation R can thus be written as: 

cos 8q — sin 6*o 



,dx H 



sin 6*o cos 9 




(3.54) 



with cos#o = X /Pi sin^o = y/p, and p = ^ x 2 + y 2 . It operates on the tensor field T 
as follows: 



(3.55) 



This tranformation law applies to partial derivatives of tensors in the same way. A 
1-dimensional interpolation will be needed to calculate fields at the points (p, 0,z) 
that are not part of the original 3-dimensional coordinate grid. This interpolation 
may require points in the x < range, where the values are calculated via the same 



coordinate transformation Eq. (3.53) as shown by the dashed lines in Fig. 3.1. 
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4.1 The TOV spacetime 

The TOV (Tolman-Oppenheimer-Volkoff) spacetime is an exact solution in general 
relativity modelling fluid balls such as isolated stars. It is the only known general 
relativistic static equilibrium configuration for spherical objects with matter. The 
line element of this spacetime is: 

ds 2 = -e u dt 2 + e x dr 2 + r 2 (d9 2 + sin 2 9dtf), (4.1) 

which contains no non-diagonal components, r in this line element represents the 
coordinate radius of the fluid ball in question, which we will henceforth call the 
Schwarzschild radius, whereas v and A are purely functions of r, indicating the time- 
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independence of this spacetime. This line element is used to model neutron stars in 
this thesis. The connection coefficients, Ricci tensor and scalar that result from this 
line element are as follows: 



GS = i«-5seH = ~ + e -»(p-7) 



Gi = Gl = Ill- \glR = \e- x (Wf - v>\' + 2 V " + 2(£ - (4.2) 

with all other components of the Einstein tensor being zero. 

Within the interior of the neutron star, the momentum-energy tensor Tj is given 

by: 

Tj = {p e +p)u l u +p5) = diag(- Pe ,p,p,p), (4.3) 

where p e here is the matter-energy density. Using GRTensorll on Mathematica, the 

Einstein field equations are computed to be as follows: 

G° - 





1 -e A 



jiii ry -ij ry 



1 x\ 

X = — + 8nre x p e (4.4) 

r 

G\: 

- 2 - e~\- + -) = -87rp 

ry £ ry £i ry 

v' = -A' + 8nre x (p e + p) (4.5) 
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r<2 _ /^3. 

Gr 2 - Cr 3 . 



I i \'\ -A/ 1 , ^ \ 1 -A // 



inp. 



(4.6) 



Eq. (4.4) can be solved for the function A when it is written as an integral as follows: 



(re" A )' = 1 - 8np e r 2 



2m(r) 



(4.7) 



where 



m(r) = Arc p e (r')r' dr' 
Jo 



(4i 



Employing Eq.s (4.4) and (4.7) in Eq. (4.5) gives us: 



. 2m + 87rr 3 p 

/ = . 

r(r — 2m) 



(4.9) 



We then use the former two equations together with the derivative of Eq. (4.5) to 



write Eq. (4.6) entirely in terms of v' and p' as follows: 



1 zA 



'-2v' + 8nre\ Pe + p))e' A (- + -) - 



2r 



[(ri/ + 1)A' + (16vrrp + 8vrrV)e A - v'\ 



inp 



P 



-{p e +p)v'. 



(4.10) 



Eq.s (4.8), (4.9) and (4.10) form a set of ordinary differential equations that can be 



solved to yield our theoretical stellar models. In our simulations, we choose a value 
for p e (r = 0) with m(r — 0) = 0, employ an equation of state giving p as a function 
of p e , and integrate the former equations from r = to a radius where the pressure 
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p vanishes. This radius defines the surface of the neutron star, which we will denote 
henceforth as R. We employ an equation of state popularly used in neutron star 
simulations called the polytropic equation of state coupled with a T law which sets 
p = np r = p e e(r — 1) for the initial data with k and T as arbitrary constants. These 
solutions give us models of the interior of the neutron stars. 

For the exterior region of the neutron star, p = and m = M where M is the 
gravitational mass of the star as evaluated at spatial infinity, which we will denote 



henceforth as the ADM (Arnowitt-Deser-Misner) mass. Therefore, Eq. (4.7) becomes: 



2M 



(4.11) 



and Eq. (4.9) yields: 



2M 



r(r-2M)' 

The line element for the exterior region thus becomes: 



(4.12) 



ds 2 



1 )dt 2 + 

r 



Y'dr 2 + r\de 2 + sin 2 <?# 2 ), (4.13) 



which is the standard Schwarzschild line element. 

Since we perform these calculations within a 3-dimensional Cartesian grid frame- 
work, we use an alternative form of the Schwarzschild line element as follows: 



ds< 



.(1 _ ™) d1 ? + (i _ ^-yV + r\dQ 2 + sin 2 a > 2 

-Adt 2 + B(dr lso 2 + r iso 2 (d9 2 + sin 2 ° ' 2 1 

-Adt 2 + B(dx 2 + dy 2 + dz 2 ), 
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where it can be computed that: 



A 
B 



t 2r iso -M 2 



y 2r 



M , A 



2r- 



(4.15) 



and: 



dr 



M 
2r- 



r iS o{2r iso + Mf 
(2rj SO ) 2 

Jr{r - 2M) 
dr; 



(4.16) 



Due to Eq. (4.11), the differential for the Schwarzschild radius in Eq. (4.16) for the 



stellar interior can be written analogously as: 



\/r(r — 2m) 
dr = — dr; 



(4.17) 



The set of ordinary differential equations for the stellar interior thus becomes: 



dp 
dr 



dm Anp e r • yj r{r — 2m) 

2{m + 47rpr 3 ) 
y/r(r - 2m)r iso 
(m + Airpr 3 ) 
yjr(r - 2m)r iso 



dv 
dr^gQ 

-(Pe+P) 



(4.18) 



The solution in the stellar interior is matched to the solution in the exterior to obtain 
a global solution. These solutions were first considered by Tolman, Oppenheimer and 
Volkoff in 1939 |80j,[66j to model isolated stars, henceforth the name TOV spacetime. 
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4.2 Perturbations on the TOV spacetime 

In this section, we present small amplitude perturbations to the TOV spacetime as 
first introduced by Thorne and Campolattaro in 1967 [79]. These perturbations are 
viewed as vector displacements of the fluid elements with respect to the coordinate 
system of the spacetime. In terms of the line element, the perturbation is introduced 
as follows: 

ds' 2 = ds 2 + h^dafdx", (4.19) 
where ds' 2 is the perturbed line element, ds 2 is the original line element of the TOV 



spacetime as described in Eq.s (4.1) and (4.13) and are the components of the 



perturbation metric. The components of the perturbation metric and the fluid ele- 
ment displacement vector can be categorized into quantities that transform as scalar 
fields, vectors or tensors under a rotational group. The quantities that transform as 
scalar fields are constants under a rotation group. They can be expanded in terms 
of scalar spherical harmonics Y^(#, <ft) and possess even parity. Those that transform 
as vectors, eg. = t^A u , are expanded in terms of vector spherical harmon- 
ics ty l m j = djY^iO^cj)) in the even parity and $^ = e^dkY^O, 0) in the odd parity 
where e\ = e| = sin# and e| = e| = 0, whereas those that transform as ten- 

sors, eg. A^ u = i^T^A mn , are expanded in terms of tensor spherical harmonics 
^mjk = Ym-jki a covariant derivative with respect to the 3- metric 7^ for the TOV 



spacetime, & mjk = j jk Y^ in the even parity and £ l mjk = \{e^ l mnk + e^¥ mnj ] 
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the odd parity. £ r , hoo, Hq t and h rr are quantities that transform as scalars; whereas 
(£0, £</>)> (hoe, ho^), (h r g, h r ^) transform as vectors, and hjk with (j, k = 9,4>) transform 
as tensors. Therefore, we can write the overall odd-parity harmonics as follows: 

( r = 0,b = U{r,t)& ma ,b = V(r,t)& n3i 

hoo = hor = h rr = 
h 0j = h {r,t)$ l mj , h xi = hi{r,t)$ l mj , h jk = h 2 {r,t)x l mjk . (4.20) 
The overall even-parity harmonics can be written as follows: 

e r = X{r,t)YHe = V{r,t)¥ m2 ,^ = V(r,t)¥ m3 

h 00 = e v Ho(r,t)Yl,h 0r = H x {r,t)Y^h rr = e v H 2 {r,t)Y l m 
h 0j = ho{r,t)¥ mj , h Xj = h!{r,t)¥ mj , h jk = r 2 K(r,t)V l mjk + r 2 G(r, t)$ l mjk . (4.21) 



Eq.s (4.20) and (4.21) can be simplified via small coordinate tranformations x M 



x v _|_ rf{x) on the perturbation metric components as follows: 

h'w = Kv + + Vvw ( 4 - 22 ) 

where r] = rj r = and rjj = A(r, t)& m , , for the odd-parity harmonics, and r] = 
M (r, t)Y^, rji = M(r,t)Y^ and rjj = M 2 (r, t)^ l m j for the even-parity harmonics. For 
the odd-parity harmonics, the simplification is performed by setting A so as to annul 
the function h 2 , whereas for the even-parity harmonics, M Q , M\ and M 2 are chosen 
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h sin 6Y l 2 



h 



(4.23) 



so as to annul h , hi and G. With these gauge simplications, the perturbation metric 
for the odd-parity harmonic with m = becomes: 



h^mOY^ 



h sin 6Y \ 2 tusmOY^ 
with the corresponding fluid element displacement vector as follows: 






C/sin0Y o ' 2 



,a = (1,2,3). 



(4.24) 



The perturbation metric for the even-parity harmonics for m = becomes: 



hfxu 



H e» 










#i 
















r 2 if 














r 2 sin 2 9K 



yl 



(4.25) 



with the corresponding fluid element displacement vector as follows: 



XY l 



02 



(4.26) 



These general forms indicate that the odd-parity harmonics are described by differen- 
tial rotations that do not change the profile of the star's density, pressure and shape, 
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whereas the even-parity harmonics are described by oscillations of the star's density 
and pressure. 

Armed with the previous general forms, we now consider the most basic pertur- 
bation mode, ie. the I = mode. The odd-parity harmonics for this mode are: 



and the even-parity harmonics are: 

£ r = XY °; V° oj = 0; *° 0jk = 



(4.27) 



(4.28) 



Eq. (4.27) indicates that the / = mode contains no odd-parity harmonics, and 
thus does not possess any differential rotations. However, the perturbation metric for 
even-parity harmonics for the I = mode becomes: 



h 




and the fluid element displacement vector becomes 



„2 e - a /2 Xy 



H e» 


H x 


Hi 


H 2 e 















(4.29) 








(4.30) 
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We can thus write the perturbed line element for the I = 0, m = mode as follows: 



ds 2 = e u (H Y° - l)dt 2 + 2H 1 Y°dtdr + e\H 2 Y° + l)dr 2 + r 2 (d6 2 + sin 2 9d0 2 ). (4.31) 



The 4- velocity of the fluid elements corresponding to the aforementioned displacement 



vector is: 



ur 



e _ v/2(1 + 

-r- 2 e-( v+x V 2 X i0 Y ° 





,P = (0,1,2,3). 



(4.32) 



From the displacement vector, we also compute the Lagrangian change in the number 
density of baryons contained in a certain stellar volume, as follows: 



An 

n 



-(& + 



2 



l e- x ' 2 X 1 Y* 



H2Y ° 



(4.33) 



where the subscripted colon indicates a covariant derivative with respect to the per- 
turbed metric to first order in the perturbation functions. The corresponding Eulerian 
changes in the density and pressure of mass-energy of the star hence become: 



An 

n 



Sp e = (p e +p)— + p e , 1 r- 2 e- x / 2 XY ° 



dp 



n 



o • 



(4.34) 



where 7 



p c +P pa 

V Pe,l 



. We note that Eq.s (4.31) to (4.34) constitute spacetime and 



fluid perturbations that do not violate the spherical symmetry of the equilibrium 
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TOV background. Substituting these expressions into Einstein's field equations and 
simplifying by neglecting higher orders of the perturbation functions, we now write 
the equations of motion for the / = mode perturbation as follows: 
5G\ = 8tt6T}: 

2 r - 1 e- {x+u) H fi - r-V A # ,i - r" 2 e~ A (l + rv')H 2 - 

87T 7 p(r- 2 e- A/2 X i - H 2 /2) - 87rr~ V A/2 p,iX = (4.35) 

SG° = 8tt<5T °: 

- r- x e- x H 2 ± + r~ 2 e- A (-l + r\ x )E 2 + 
87r(p e + p)(r- 2 e- x/2 X A - H 2 /2) + 87rr-V A/2 p e>1 X = (4.36) 

SG 2 = 8tt5T 2 , 5G 3 3 = 8tt5T 3 3 :: 

- e-"^ 2)0 o/2 + e'^H lfil + r - x e~ x+v \2 - r\ tl )H lfi /2 - 
e- x H 0A1 /2 - r~ V A (2 - r\ A + 2rz/ 1 )# , 1 /4 - r - 1 e~ x (2 + ru^H^/4 + 

4r- 1 e- A ((A, 1 - i/ )1 )(2 - r Kl ) - 2rv M )H 2 - 
8vr 7 p(r- 2 e- A/2 X i - H 2 /2) - 8nr~ 2 e- x/2 p A X = (4.37) 

SGI = 8tt 5T^. 

r- l e~ x H 2fl - 87r(p e + p)r-V A / 2 X = (4.38) 

5G\ = 8ir6T?: 

r-'e-^^A.i + i/,i)ifi - r^e- v H 2fi + 8vr(p e + p)r-V" +A/2 X = (4.39) 
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The Euler equation for this perturbation mode, which is the projection of 5T±. U = 
along the direction orthogonal to the 4-velocity of the fluid u, is as follows: 

(Pe + p)Cooe A+ ^ - { -^^H 0A Y ° + (S Pe + Sp) ^ = -Sp,i. (4.40) 



Eq.s (4.35) to (4.40) can be solved to yield a first order ordinary differential equation 



for the perturbation function X. To do this, we write the perturbation function in 
terms of a radial dependence and a harmonic time dependence, X(r,t) = X(r)e %UJt . 



We then substitute Eq.s (4.34), (4.35) and (4.38) into (4.40) to obtain the following: 



)(r- 2 e- A / 2 X n + - Z + A7r nP e x - ^)X A + 

Pe,l Pe,l 2 



+ — e A - Z A - 4vr(p e + p)Zre x + u 2 e x ~ v )X = 



(4.4i; 



We now consider the I = 1 mode perturbations. The odd-parity harmonics for 



this mode are: 



$02 = 0;$i 12 



3 l '\ 



±id 



$03 



^-) 1 / 2 sin 2 0;<l>i 1 3 = T^^V^smflcosfl 
Xmjk = 0' h k = 2,3, 



(4.42) 



whereas the even-parity harmonics are: 



1 / ^ \1/ 



o-V^ /2 cos6;Yl 1 = TK - 



3 l '\ 



e^sinfl 



* 02 = "(f ) 1/2 sin^i 12 = T(h 1/2 e^cos6 
^03 = 0;*I 1 3 = -^(^) 1/2 e ± ^sin^ 
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$m22 



.vi/ 1 — yi. (j) 1 

^m22 _I mi v m33 



-^33 = sin 2 fl^; $ m ^ 3 = = 0. (4.43) 



In the spirit of Eq. (4.20) and performing a simplification via small coordinate trans- 



formations previously done for the / = mode, which sets G — K = 0, we can then 
write the metric for the odd-parity I = 1 perturbation mode as follows: 



h 



H e u iuHi 
iuHt H 2 e x 











Y 1 

mi 



(4.44) 





where we include a harmonic term in the 10— and 01— components in order to 
simplify the resulting system of ordinary differential equations shown further on. 
Correspondingly, the fluid element displacement vector is: 



r 



-2 -A/2 W 1 



r 2 m 



(4.45) 



where we see that W is the perturbation function describing radial fluid oscillations 
whereas V describes azimuthal fluid displacements. The perturbed line element for 
this perturbation mode is thus: 

ds 2 = e u (H Y^-l)dt 2 + 2iuH 1 Y^dtdr + e x {H 2 Y^ + l)dr 2 +r 2 {d6 2 + sm 2 6d<j) 2 ). (4.46) 

With this line element, the Lagrangian change in the baryonic number density be- 



comes: 



^ = (rt-^Xi + 2r" 2 V - ^)Yl 



(4.47) 
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which results in the Eulerian changes in the pressure and matter-energy density as 
follows: 

5p e = {^±2l e-x/'Srt + p^r-'e-WWY* 

5 P = e-^SiY* + ir-V A / 2 WY^ (4.48) 

where Si = ^pe v / 2 {r~ 2 e~ x l 2 W t i + 2r~ 2 V — ^). To first order in the perturbation func- 
tions, again using GRTensorll on Mathematica, we obtain the Einstein field equations 
as follows: 
5G\ = 8tt5T}: 

(2r- 1 e- x -"H lfl - r^e^H^ - r~ 2 e~ x H 2 (l + ru A ) + H r~ 2 )Y^ - 

87re~ u/2 Si - %np,ir- 2 e- x/2 WY^ = (4.49) 



SG 2 = 8tt5T 2 : 

[ g ( " 2 + rX ^ Hl >° ~ ~Y H ^o + —j— (-2 + rA,i - 2rui)H 0> i 



(2 + rui)H 2 ,i + e- x -"H lfil - ^#o,n + 
^(e-V((A i - K i)(2 + rui) - 2ri/ u ) - 2) + ^]^ t - 



STre^Si - 87rp 5 ir- 2 e- A/2 WY* = (4.50) 



<JG§ = 8tt<5T : 



[-r^e-^,! + r " 2e 2 " Ag2 (2(-l + r\,i) - 2e x )\Y^ + 



s J Pe+P K -"/ 2 SiY^ + 87r Peil r- 2 e- A / 2 Wj, = (4.51) 



IP 
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5Gl = 8tt6T^. 



x eT x E 2 ,i + r~ 2 e~ x H 1 - 8n(p e + p)r- 2 e~ x/2 W fl = (4.52) 



SG J = 8 ^ 5T o : 



r 2 , „ i „ e A 



(# 2>Q - e- x H x ,x + — (A,i - v tl )H x - 87r(p e + p)r"V = (4.53) 
z, z 

8G{ = 8n5T(: 

H lj0 = e v H Q>1 - e\r- x - ^)H + e u {r~ l + ^)H 2 . (4.54) 



We write the perturbation functions in terms of a radial dependence and a harmonic 
time dependence, ie. H (r,t) = H (r)e iuJt , H x {r,t) = iJ 1 (r)e^*, H 2 (r,t) = H 2 {r)e iuJ \ 



V(r,t) = V(r)e iut and W(r,t) = W(r)e iujt . We then decouple Eq. ( |4.51[ ) into two 
equations for Hi and H 2 as follows: 

H 2 = -r- 1 ^ + r -1 87r(p e + p)e A/2 iy (4.55) 

= ^(^(A i e A )#i + r _1 87r(p e + p)e 3A/2 iy - 16vr(p e + p)eV. (4.56) 



Eq. (4.54) however can be simplified to yield: 



rH 0>1 = -(— i - 1)# - (^r 1 + 1)#2 + iojre~ u H 2 . (4.57) 



We also consider the perturbation of the energy-momentum conservation as follows: 

KT^) = o 

- (p e + p)e~ u u 2 V = -p A r- 2 e- x/2 W - mr~ 2 e- x/2 W yl + 

fH 2 + ^H -2 1P r- 2 V, (4.58) 
62 



Chapter 4 



Stellar perturbations and phase transitions 



which can be decoupled into two equations for V and W as follows: 



167rrV(p e + p)e x ~ v V = -87r(p e + p)i/ tl e v/2 W + - 2ru; V)iJi + 



(167rr 2 p e e A - 3r\ yl )H 



(4.59) 



Swpje^W^ = -ru 2 e" , H 1 - (^ - 4n P1 r 2 e x )H 2 - 



(1-e 



#o - 87rp il e A/2 W r - 167cjrye x V. (4.60) 



Eq.s (4.55), (4.56), (4.57), (4.59) and (4.60) form a system of third order ordinary 



differential equations for the even parity / = 1 perturbation mode and can be solved 
to obtain the mode frequency u. To consider the effect of this perturbation on the 



exterior TOV background spacetime, we solve Eq.s (4.51), (4.54) and (4.58) in the 



vacuum surrounding the star. Using Eq.s (4.11) to (4.13), we obtain, as in [19J: 



1 s 3 /3 + 8M 2 (3 



oo 



3 <r(l-<r) 



2Mq 



1-0 



rAo! H 2 



<0 2 A 



(4.6i; 



where q = — and (3 =arbitrary function of the coordinate time. Considering in- 



finitesimal coordinate transformations, we use Eq. (4.22) to obtain the following for 
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the even-parity I = 1 mode: 

hw = [e u H -2M 0fl + u 1 e u - x M 1 }Yi 

hw = [^-Mo^+utMo-M^Yl 

h vv = [e»H 2 -2M 1A + \ A M 1 }Yi 

h , f = -[M 2 , + M )¥ mj 

h Vj , = -[M 2>1 -2r- 1 M 2 + M 1 ]¥ mj 

h rk , = -2[re~ x M l + M 2 ]¥ mjk . (4.62) 



Eq.s (4.62) show that for the even-parity / = 1 mode perturbation metric to preserve 



the same form, the following will have to hold true: 

ho'j/ = hi'ji = hjiy = 0, (4.63) 

which entails: 

M 2j0 + M = 0, M 2 ,i - 2r~ 1 M 2 + M 1 = 0, 2re~ x M 1 - 2M 2 = 0. (4.64) 



Eq.s (4.64) can be solved to yield the following: 

M = -a f, M 1 = ar~ 1 e x f 1 M 2 = a/, (4.65) 
with / = r exp[J r °° r _1 (l — e x )dr] and a =arbitrary function of the coordinate time. 
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Substituting Eq. (4.65) back into the first three equations in (4.62), we obtain: 



H + (2a fi0 e- u + ar~V)/ 



H[ = H l + a [2 r -\l-e x )-v l ]f 



H> 2 



W 
V 
S' 



J ff 2 -ar- 1 [2r- 1 (l-e A ) + A, 1 ]/ 
W - are x/2 f 
V + af 
S. 



(4.66) 



Therefore, the perturbation functions H , Hi, H 2 , V and W are only unique up to 



the transformations shown in Eq. (4.66). Given Eq.s (4.66), the perturbations in the 



exterior background spacetime Eq.s (4.61 ) can be set to zero which then preserves its 



spherical symmetry. When H' Q = 0, we obtain a(t) = —^M/3t. Using this gauge not 
only guarantees that the exterior background spacetime is spherically symmetric but 
as a unique gauge, it also can be used to remove the gauge arbitrariness in the former 
perturbation functions. 



4.3 Non-radiative pulsations of neutron stars and 
their frequency modes 

In the previous section, we have delved into non-radiative perturbations on the TOV 
spacetime, ie. the / = and even parity / = 1 modes. In this section, we will apply 
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the formalism in our neutron star models to obtain their non-radial pulsation modes. 
By considering boundary conditions at the center and at the surface of the star, we 
employ a shooting method to solve the systems of ordinary differential equations in 
both the I = and even parity I = 1 modes obtained in the previous section. 
For the I = mode, we follow Misner et al |60j as well as Kokkotas and Ruoff 



in re-writing Eq. (4.41) as: 



where: 



±(pf) + (Q + u J 2 W K = 0, (4.67) 
dr dr 



( = r 2 e~ v X, (4.68) 
r 2 W = (p e + p)e^ 3X+ ^ 2 
r 2 P = 7 pe (A+3 " )/2 

r 2 Q = e (A+3,)/ 2(pe + + 2^1 - 8 7re 2A p). (4.69) 

4 r 



Eq. (4.67) can itself be decoupled into a system of two ordinary differential equations 



as follows: 

dr P 

^L = -(u*W + Q)C. (4-70) 
dr 

We now consider the boundary conditions for neutron star models. At the center of 
the star, the radial displacement of the fluid element vanishes, whilst at the surface 
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of the star, the Eulerian change in the pressure of the star vanishes. These conditions 



translate respectively as: 



C(r = 0) = 0, 



(4.71) 



TUT E)3 

(r- 2 e- x/2 W) tl = ( 7j R)- 1 (4 + {-^)e x W + u 2 {^-)e- u )(r- 2 e - x / 2 W), (4.72) 

R M 

where R denotes the radius and M the ADM mass of the star [9]. [49] found that 
via a Taylor expansion, £(r) ~ ($r 3 + 0(r 5 ), and r)(r) ~ r] + 0(r 2 ), hence we obtain 
C,i ~ and rj^i ~ as r — > 0. We then employ a simple shoot and match method 



to solve the boundary value problem Eq.s (4.70) to (4.72) to obtain the normal 



modes of pulsation for the neutron star. We input a range of test u 2 values and for 



each of these values integrate Eq. (4.70) from r = to r = R using a simple finite- 



differencing scheme and observe the difference between the integrated value and the 



value imposed by the boundary condition Eq. (4.72) at r = R. The u> values that 



give us zero difference are the normal mode frequencies for the neutron star. We note 
that using this method, convergence is achieved very rapidly. We denote the lowest 
of these frequencies as the fundamental mode. 

For the even-parity / = 1 mode, we similarly consider the boundary conditions at 



both the center and at the surface of the star for Eq.s (4.56), (4.57) and (4.60). We 



note that at the center of the star, the perturbation functions Ho, Hi and W must be 
finite whilst at the surface of the star, they must match with the values that preserve 



the spherical symmetry of the exterior background spacetime. These conditions thus 
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translate as follows: 



H ~ hr; Hi ~ —8ii(p e c + p c )wr 2 ; W ~ wr 2 



(4.73) 



where ft, and w are arbitrary constants, and p e c and p c are respectively the matter- 
energy density and pressure at the center of the star, and: 



H (R) =0;H 1 {R) = 0. 



(4.74) 



This boundary value problem Eq.s (4.56) to (4.60), and Eq.s (4.73) to (4.74) can 



then be similarly solved using a shoot and match method where test u 2 values are 



chosen and Eq.s (4.56) to (4.60) are integrated from the center to the surface of the 



star until a match occurs with the boundary conditions at the surface of the star. In 
choosing the test u 2 values, we follow [56] in using a variational principle for the I = 1 
mode frequencies as obtained by Detweiler in 1975 



LO' 



[ e -^ 2 [p e +p)e 
Jo 



A 2 £ + 2 l/ 2 )-e- A/2 P* 
r 2 87r 



Pe+P 167T 

p e RMW 2 {R) 
R 1 ' 



(4.75) 



We employ test functions Hi = — 87r(p ec + p c )r 2 [l — (r/R) 2 } and W = r 2 [l — (2r/i?) c 



and substitute Eq. (4.47), Eq. (4.55) and (4.56) into Eq. (4.75) to obtain the test 



bj 2 values for the neutron star. We note that the variational principle approach is able 
to give us a good estimate of the mode frequencies even without an exact knowledge 
of the perturbation functions. 
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For these perturbation modes, when the fundamental u 2 obtained is negative, the 
star's oscillation increases without bound and the star is said to be unstable. When 
the fundamental mode is zero, the star is at a critical point between a stable branch 
and unstable branch. The star will remain at this point unless perturbed. We shall 
delve into this in the following section. 

4.4 Equation of state change of neutron stars and 
its instability time scale 

In the previous sections, we present the formalism and the methods by which we 
obtain the time scales associated with non-radiative pulsation modes for non-rotating 
neutron stars. In this section, we consider the time scales involved when a non- 
rotating neutron star undergoes phase transitions induced by a change in its equation 
of state. Since an analytic result has not been established governing how the time 
scale of collapse of a neutron star varies with the speed of its phase transition, in this 
section, we shall consider numerical results. 

We consider isolated static neutron star models with a polytropic equation of 
state p = Kp r as mentioned in Section 4.1, with k = 80. The rest mass of these 
neutron stars are set to be at the maximum allowed for a static equilibrium star 
with the adiabatic index T = 1.9. According to Harrison et al [41] . the adiabatic 
index has the upper limit of 2 imposed by the causality condition. It is a well-known 
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Figure 4.1: 



Newtonian result that the speed of sound has the ultrarelativistic limit f5 soun d = 
(jfc)V2 = [ e ( r _ i)]i/2 (r - l) 1 / 2 . Therefore, we set T = 1.895 in order to allow 
a change of equation of state with T increasing but not exceeding the upper limit 



of 2. Fig. 4.1 shows the variation of the rest mass and ADM mass with respect to 
the central matter density of stars modelled with equations of state with different 
adiabatic indices. As mentioned in the previous section, there is a maximum for each 
of the curves, which indicates a critical point between the stable branch on the left 
and the unstable branch on the right. We see that as V increases, the curve moves 
further to the right with decreasing maximum. As the central matter density is less 
than 1 in the geometric units used, the pressure of the star as well as the number of 
baryons packed in the star decreases as T increases. 
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Figure 4.2: 



We set the rest mass of the stars to be equal to the maximum rest mass of a star 
with T = 1.9, ie. at Mf, = 2.2135. With this setup, the stars are thus set on the 
stable branch of the T = 1.895 curve, with their phases allowed to transition with 
time to the maximum point of the T = 1.9 curve. We utilize the GRAstro-2D code 
and use a finite-differencing resolution of dx = 0.025. Convergence of the results with 
respect to resolution is verified by performing simulations using other resolutions, 
eg. dx = 0.05, 0.075. A slow change is imposed on the equation of state by slowly 
increasing T at a constant rate simultaneously conserving baryonic mass, which 



is checked to be conserved to the level 0.00074%, shown in Fig. 4.2 Two categories 



of cases are considered, namely: 



(a)r = r + — t,t < t n (b)T = r + -gb,t < t f2 . 
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In (a), T increases to T — 1.9, reaching the critical point of the rest mass-matter 
density curve for V — 1.9 at t — t/i, thus becoming unstable for t > tf\. In (b), 
T increases to Y = 1.905, overshooting the critical point beyond which there is no 
equilibrium configuration with the given rest mass, for V > 1.9. Within these two 
categories, neutron stars undergoing three rates of change of are investigated, 
namely: 

<i) = 0.00005, (ii)— = 0.000075, (in)— = 0.0001. 
at at at 

We determine the time scale of the collapse of these neutron stars as they cross 
the critical point by measuring 5p = p — pi, where i is the central matter density of 
the star at the critical point of the rest mass-matter density equilibrium curve for 
T = 1.9, and by measuring the rate of change of dp and its acceleration. 



Fig. 4.3 shows the change of Sp, dSp/dt and d 2 5p/dt 2 with time for each of the 
transition speeds (i),(ii) and (iii). From this figure, we observe that 5p changes at a 
constant speed for each of the transition speeds, with d 2 5p/dt 2 ~ 0. 



Fig.s 4.4 to 4.6 show the comparison of the Sp, dSp/dt and a I Sp/dt changes with 
time, between scenario (a) and (b). From these figures, we observe that the changes 
of 5p and their respective rates of change do not differ much between a neutron star 
that undergoes a phase transition up till the threshold point, and a neutron star that 
undergoes a phase transition past the threshold point reaching a state which cannot 
be described by an equilibrium configuration. Whilst a small difference is seen only 
toward the late stage of the collapse, the difference in the d 2 5p/dt 2 is negligible when 
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Figure 4.5: 



74 



Chapter 4 



Stellar perturbations and phase transitions 




75 



Chapter 5 

Critical gravitational collapse 

5.1 Critical phenomena in general relativity 

Critical phenomena in general relativity were discovered by Choptuik (20] in 1993 in 
numerical simulations of spherical scalar field collapses. The phenomena of universal- 
ity, mass scaling and self-similarity observed in these gravitational collapses garnered 
the name of critical phenomena, analogous to phase transition phenomena observed 
in condensed matter systems. 

In particular, as a parameter of a set of general relativistic initial data is varied, the 
evolution of this initial data passes through a threshold between black hole formation 
and dispersion to infinity. They evolve toward a spacetime which can be stationary 
or scale-free [38]. We shall call this spacetime the critical set henceforth. These 
evolutions leave the threshold via the one unstable mode of the critical set, and the 
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time scale of this departure is denoted as the critical index. In the present context, 
universality means that the critical index is independent of the initial data parameter 
that is varied. 

Mass scaling occurs when the black holes that form at this brink of collapse begin 
with infinitesimal masses that scale with respect to the distance of the initial data 
from the threshold, as follows: 

Moc(p-p*y, (5.1) 

where M denotes the mass of the black hole, p denotes the parameter of the initial 
data that is varied, p* the value of this parameter when the initial data evolution 
stays on the threshold as time goes to infinity, and 7 the critical index. According 
to [38], an initial data set in general relativity can be denoted as a function of the 
spatial coordinates, eg. S(x) where we express it in a 1-dimensional spatial coordinate 
system, and its evolution in the coordinate time can be denoted as S(x, t). S(x) can be 
the density distribution of the matter configuration or the 3-metric of the spacetime, 
as mentioned in Section 3.4. Denoting the critical set itself as S*(x), a solution of the 
initial data S(x) that produces an evolution very near the critical set can be linearized 

in a small neighborhood of the critical set as follows: 

00 

S(x,t) « S*(x) + Y,C k ( P )e x * t S k (x), (5.2) 

fc=0 

where the A^'s are eigenvalues of the critical set which can be purely real, complex 
or purely imaginary. As the critical set by definition possesses one unstable mode, as 
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t — > oo, all perturbations about the critical set vanish leaving one with positive real 
Afc, for which we denote the amplitude of the constant as Co(p*). We can then further 
linearize in a small neighborhood about the critical set via a Taylor's expansion as 
follows: 

lim S(x,t) « S*(x) + —^(p-p*)e Xot So(x). (5.3) 
t^oo dp 



We can choose a t = t* where Eq. (5.3) still holds and denote 



S(x,t*) &S m (x)+eS (x), (5.4) 



where: 

e =^( p _ p *) e Aot* (55) 

dp 

A black hole that forms in the solution S(x,t*) starts with a mass that grows 



proportionately to e **, which, in line with Eq. (5.5), is in turn proportional to 
(p — p*y/ x ° = (p — p*) 1 . From Eq. (5.5) too, we can calculate the critical index as: 



1 

A in |p — p * | 

where we have taken C = ln(e/^)/Ao. Analogous to critical phase transitions in 
materials, ie. first and second phase transitions, with discontinuous and continuous 
order parameters, we can categorize critical gravitational collapse scenarios between 
those that exhibit mass scaling and those that do not. Scenarios that do exhibit 
mass scaling are called Type II scenarios and those that do not are called Type I. 
Type II scenarios occur when the system in question does not have a preferred scale, 
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transformation 
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Figure 5.1: 

eg. in scalar field systems, whereas in Type I scenarios, there exists a scale that 
typifies the Einstein field equations for the system in question, and which cannot be 
neglected. Therefore, in Type I scenarios, the mass of the black holes formed start 
with finite value. The critical index in Type I collapse scenarios is similarly defined 



as in Eq. (5.6). 



The next feature seen in critical gravitational collapses reminiscent of phenomena 
seen in condensed matter systems is the feature of self-similarity. Self-similarity de- 
scribes the symmetries of the critical set of the collapse scenario, where the system 
exhibits scale-free behavior. There are two types of self-similarity observed in collapse 
scenarios, namely discrete self-similarity (DSS) and continuous self-similarity (CSS). 
Cahill and Taub [18] first computed the equations for CSS of a spherically symmetric 
configuration of perfect fluid, as follows: 



^9a/3 — 2# 
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where £ is a homothetic Killing vector producing this self-similarity transformation, 
which in terms of the coordinate system x^ 1 = (t, x l ), can be denoted as: 



d_ 
dt' 



(5. 



and g a s is the 4-metric of the spacetime exhibiting CSS. Eq. (5.7) reduces Einstein's 



field equations to a set of ordinary differential equations. The resulting equations for 
the spacetime in such a collapse scenario are: 



^-(.^uM ~ A? — 0) ^Gap — 0. 



(5.9) 



This leads to the following equation for the energy-momentum tensor: 



£{r a p — o, 



(5.10) 



which, coupled with Eq. (5.7), produces, in the case of a perfect fluid, the following 
for the matter variables: 



C{P = —2p; C^U a = -U a ] C^Pe = ~2p e . 



(5.11) 



Eq. (5.9) describe geometric self-similarity whereas Eq. (5.11) describe physical self- 



similarity, which is equivalent to geometric self-similarity only in the case of perfect 



fluids. Eq.s (5.7) and (5.11) can be viewed as a similarity transformation in general 



relativity whereby a spatial hypersurface is pushed forward into itself as depicted in 



Fig. 5.1| (refer to Appendix B). DSS is a more general form of self-similarity [2S] and 
the 4-metric g a p for the spacetime undergoing such a collapse scenario can be written 
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in the form: 



(5.12) 



Due to the scale-invariance of Type II critical gravitational collapses, the self-similarity 
feature can only be observed in Type II scenarios. 

Since the first discovery of critical gravitational collapse by Choptuik in spher- 
ical scalar fields, universality, mass scaling and self-similarity have since been ob- 
served in a variety of other systems such as perfect fluids [29], [62], 2D sigma models 
P|.|53j.|5ij. SU(2) models [22] , [12] , [13] , [23] , [59] and primordial density fluctuations 
in cosmological models [S3] , [SS] , [S2] , [3E] , [32] • However, these systems are largely re- 
stricted to spherical symmetry and require fine tuning of initial data and are therefore 
less prevalent in realistic astrophysical systems. Later on in 2003, Choptuik and his 
collaborators ventured into the study of scalar fields in axisymmetry [24|. Abrahams 
and Evans studied axisymmetric collapse in vacuum gravity [Tj and Lai [51] found 
Type I critical phenomena in boson star systems. Noble and Choptuik then reported 
Type II critical phenomena in spherically symmetric static neutron stars [63J. In 
2007, Jin and Suen [17] produced evidence that Type I critical collapse occurs in 
axisymmetric neutron star systems. A crucial finding in the latter that differs from 
the previous ones is that when the adiabatic index of the equation of state is varied, 
the same Type I critical collapse is seen to occur in the neutron star system. This 
presents evidence that critical phenomena not only occur via fine-tuning of initial 
data, but can occur as the equation of state of realistic star systems softens and when 
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this time scale of equation of state change is longer than the time scales involved in 
the critical collapse. 



5.2 Universality and the critical index 

In the previous section, we reviewed the basic concepts of critical gravitational collapse 
and the gravitational systems in which the phenomena has been observed. In this 
section, we focus on several aspects of the concept of universality, which drive the 
main analyses work in this thesis. 

As mentioned in the previous section, critical collapses are said to be universal 
as they produce the same critical index even when different parameters of the initial 



data are varied. This is seen to be the case due to Eq. (5.3) where only one growing 
mode survives in the critical solution carrying the system away from the critical set 
- the main criterion for a solution to be critical. Let us take a certain set of initial 
data that is tuned via variation of a certain parameter to form a critical solution. We 
then fix this parameter and vary another parameter to tune the initial data until it 
produces a critical solution. We note that these two initial data sets that produce 
the critical solution lie very close to each other. Therefore, even when the critical 
index extracted from these two sets of 1-parameter variations are the same, we cannot 
conclude that we observe universality in the strictest sense. 

Another aspect of universality is the characteristic of the critical solution whereby 
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all details of the initial data, except for the distance of the initial data from the 
threshold, and its conserved geometric quantities, are washed away in the dynamical 
evolution, thus arriving at a universal configuration. In scale-invariant spacetimes 
where self-similarity is observed, universality in this aspect indicates that the sys- 
tem has entered a region of intermediate asymptotics which simultaneously ceases to 
depend on details of the initial data as well as being far away from equilibrium. In 
Type I collapses of neutron star systems where there exists a constant entropy in the 
evolution, it is especially unclear how dynamical processes drive the system toward a 
universal configuration [38J. In light of this, observing evolutions of initial data sets 
that are very close to each other does not enable us to strictly conclude universality. 

5.3 The dynamical systems picture 

Another approach in understanding critical gravitational collapse which proved to be 
very useful, is the dynamical systems picture. The theory of dynamical systems is 
generally employed in studying physical systems that evolve in time, where the state 
of the system at an instant in time is described by an element x in a phase space X, 
which can be finite or infinite-dimensional. The evolution of the system is described 
by an autonomous system of differential equations on X as follows: 




(5.13) 



S3 



Chapter 5 



Critical gravitational collapse 



and: 



(5.14) 



where we see that the right hand sides do not have an explicit time dependence. When 



X is finite-dimensional, Eq. (5.13) becomes an autonomous system of ordinary differ- 



ential equations. When X is infinite-dimensional, Eq. (5.13) becomes an autonomous 



system of partial differential equations. Eq. (5.14) is a smooth function which gener- 



ates a flow <f>(t, x). Given an initial condition x(0) = xq, a solution <p(t, xq) is obtained 



and is denoted as the trajectory of Eq. (5.13) based at xq. An important feature of 



this equation is the invariance of / up to translations in time, which therefore en- 
ables the translation of solutions based at to ^ to that based at to — 0. Among 
the solutions to this system of differential equations, there exists an important class 
called fixed points, where the f(x) = 0. A fixed point x is said to be stable when 
every solution <j>(t, xo) with x lying in a neighborhood of x approaches and remains 
close to x as t — > oo. The fixed point is further said to be asymptotically stable when 
(f)(t, Xq) — > x as t — > oo. Conversely, the fixed point is said to be unstable when every 



solution 0(t, xo) leaves x as t — > oo. Linearization of Eq. (5.13) in a small neighbor- 



hood of x enables one to cast the system of equations into a characteristic form and 
obtain its eigenvalues and their corresponding eigenvectors, as follows: 



dx 
~dt 



x 



x j (t) 



Ax 



e 3 v J 



(5.15) 
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where Xj is the eigenvalue associated with the eigenvector vK Solutions that lie in 
the linear subspaces spanned by the eigenvectors form what is called invariant sets. 
These invariant sets can be further and correspondingly categorized as stable and 
unstable depending on the directions of the trajectories. The eigenvalues are crucial 
in determining the nature of the fixed point as well as the trajectories in its neigh- 
borhood. If there exists more than one eigenvalue, and if all of them possess negative 
real parts, the fixed point is determined as asymptotically stable. If one of the eigen- 
values have a positive real part, the fixed point is unstable. If all the eigenvalues 
have non-zero real parts, the fixed point is called hyperbolic. Eigenvalues that are all 
purely imaginary indicate that the fixed point is stable but not asymptotically stable. 
Different combinations of the eigenvalues also produce different trajectorial behaviors 
near the fixed point. We present here a 3-dimensional example where the eigenvalues 
consist of a conjugate pair of complex numbers, and a real number. Let us take the 
system of ordinary differential equations as follows: 



x 
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-P 
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a 











A 



X 



;u — a + ip, 



(5.16) 



which has the following solution: 







X 




y 


(*) = 


z 





e at ((cos/3t)x(0) - (sin/ft)y(0)) 
e at {{sm(3t)x{0) + (cos/9t)y(0)) 
e xt z(0) 
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Figure 5.2: 



The fixed point for Eq. (5.16) is the point (0,0,0). In this case, the trajectories of 



solutions in the small neighborhood about this fixed point form a spiral structure. As 



the radius yi 2 + y 2 of the trajectory shrinks, the z component of it increases. 

Using the dynamical systems approach, general relativistic initial data sets that 
evolve and form a critical solution can be seen as points in an infinite-dimensional 
phase space. As shown in Section 5.1, linearization can be performed in a small neigh- 
borhood of the critical set, and its eigenvalues determine the nature of the trajectories 
of these initial data sets in the neighborhood of the critical set. The critical set is a 
semi-attractor that serves as a hyperbolic limit set possessing infinite decaying modes 
but one growing unstable mode in its neighborhood. Thus, trajectories are attracted 
to the critical set via decaying modes that are tangent to the trajectories in the phase 
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Figure 5.3: 



space, but as they reach the critical set, they are repelled from the threshold via the 



one growing unstable mode (Fig. 5.2). The subspaces that are inhabited by these 
trajectories are manifolds. Manifolds that contain trajectories leaving the critical 
set are unstable while those that contain trajectories attracted to the critical set are 
stable. These manifolds are tangent to their corresponding linear subspaces which 
are obtained in the former linearization in the small neighborhood of the critical set, 
and are of the same dimension as the latter [37]. The nonlinearity of the system 
is what allows for the possiblity of the critical set being a closed or periodic orbit 
rather than just a point. In the case of a periodic orbit, the solution takes the form 
(j)(t,Xo) = (p(t + T,xo) as t — > oo. As mentioned in Section 5.1, the critical set in grav- 
itational collapse scenarios can thus be either a limit cycle or a fixed point. In Type 
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II collapses, the critical set is a limit cycle when the the system exhibits CSS, and a 
fixed point when it exhibits DSS. The stationarity of the critical set thus precludes 
any possibility of its emission of gravitational radiation by definition. 

Recalling Section 5.1, in critical gravitational collapses, 1-parameter variations of 
the initial data produce groups of initial data that arrive at different end states. These 
end states are separated by a threhold into two basins, ie. that of a black hole and 
that of dispersion into infinity or a stable compact object. The closer the initial data 
point is to the threshold, the longer its trajectory will stay parallel to the threshold 
and the easier it gets attracted to the critical set. As these initial data points are only 
very near but not directly on the threshold, their trajectories are repelled by the one 
unstable mode of the critical set and leave the threshold in a path perpendicular to the 



threshold even before it reaches the critical set (Fig. 5.2). The manifold that contains 
only trajectories leaving the critical set is called the unstable manifold of the critical 
set. Points that lie directly on the threshold by definition produce trajectories that 



never leave the threshold (Fig. 5.3). This is due to the fact that these trajectories 
will all lie within the stable manifold of the critical set. We thus note that the 
critical set defines an attraction basin, a region of the phase space near the threshold 
inhabited by initial data points that produce trajectories that are attracted to the 
critical set. The size of this basin can be described in terms of its width along 
the direction perpendicular to the threshold, and its breadth along the tangential 
direction. The threshold, on the other hand, is a plane of codimension one in the 
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phase space separating between the black hole and dispersion/stable compact object 
attraction basins. We denote the region in the phase space situated between the stable 
and unstable manifolds as a hyperbolic region where trajectories are both moving 
toward the critical set and also moving away from the threshold. Strictly speaking, 
in the theory of dynamical systems, the trajectories in this region are sequences of 
points rather than a curve in the phase space. 

The critical set itself however is denoted as the center manifold, as all its eigenval- 
ues have zero real parts, in contrast with the unstable manifold in its neighborhood 
whose eigenvalues have positive real parts and the stable manifold where the eigen- 
values have negative real parts. In dynamical systems theory, the center manifold is 
where bifurcation occurs. Bifurcation describes the change in the qualitative structure 
of the solutions when parameters in the defining differential equations are varied and 
when they enter a range of values, called bifurcation values. The bifurcation point is 
the point where the qualitative behavior switches. The center manifold thus contains 
the bifurcation point, and is tangent to the center linear subspace. Codimension one 
bifurcations depend only on one parameter. We present here an example where a 
particular type of bifurcation occurs, called the Hopf bifurcation. Let us take the 
system of ordinary differential equations as follows: 



x 



y + x(l-x 2 - y 2 )(A 



2 2\ 

x -y ) 



y = x + y(l-x 2 -y 2 )(A 



2 2\ 

x -y ), 



(5.18) 
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which can be rewritten as follows: 



2 2,2 

r = x + y 



rr = xx + yy 
r = r(l-r 2 )(4-r 2 ). (5.19) 



The trajectories of solutions to Eq. (5.18) will be governed by the signs of r in 



Eq. (5.19) for different ranges or r. If the initial condition ro is such that < ro < 2, 
the solution r(t,r ) approaches the value of 1 as t — > oo. However, when r > 2, 
the solution goes to infinity as i oo. Thus, the limit cycle where r = 0, r = 2 
is a repelling set, with the direction of the trajectories changing with respect to the 
variation of the parameter r in the initial condition. We note however, in numerical 
simulations of nonlinear infinite-dimensional systems where no exact solutions exist, 
bifurcation is never actually seen, due to the fact that the evolutions cannot in princi- 
ple be carried out for t — > oo. In effect, what we see is only a jumping back and forth 
of initial data points in a small neighborhood about the threshold. This is analogous 



to the introduction of a small number |e| in Eq. (5.18). The trajectories are therefore 
seen to leave the threshold in either direction, never actually dwelling on the critical 
solution as t — > oo. 
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6.1 Time scale comparisons with non-radiative pul- 
sation modes of neutron stars 

In this section, we shall seek to clarify the different time scales involved in axisymmet- 
ric neutron star critical collapse. Considering the results obtained by Jin and Suen 
[17], we obtain the time scale of attraction of the system toward the critical solution 



and the dominant frequencies of the oscillations of the critical solution. Fig. 6A shows 
the evolution of the lapse, density and 4-scalar curvature at the center of collision. 
Via a Fourier analysis, the dominant frequencies of the oscillations of these three evo- 



lution variables are obtained. Fig. |6.2| indicates two dominant frequencies, namely, 
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Ui = 0.15 and u 2 = 0.23, corresponding to the oscillation periods of Ti = 0.21 and 
T 2 = 0.13 respectively. These values are in units of solar masses and all values hence- 
forth in this chapter will be designated in these units. These frequencies are found to 
be similar across all three variables. Using the dominant frequencies, we then fit the 
critical solution with the following equations: 
lapse: 

a = + ae~ A( *~ 90) + fc(cos(u;i(£ - 90) + c) + ratio a cos(w 2 (t - 90) + d)) (6.1) 
density: 

p = l — , ) - 6(cos(wi(t - 60) + c) + ratio p cos(w 2 (t - 60) + d)) (6.2) 



4-scalar curvature 
density x 50 
lapse 




50 100 150 200 250 300 

t 



Figure 6.1: 
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Figure 6.2: 
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4-scalar curvature: 

R = 1 im J - 6(cos(a;i(t - 100) + c) + ratio R cos(u 2 (t - 100) + d)), (6.3) 

— h ae~ Alt ~ luu - ) 

where a, b, c and d are fit constants, aoo, and Roo are respectively the averages 
of the lapse, density and 4-scalar curvature throughout the oscillation phase of the 
critical solution, ratio a = 0.2, ratio p = ration = 0.4 are the ratio of the amplitudes of 
u>2 over u>i obtained from the Fourier transforms of the respective evolution variables. 



The fit results are shown in Fig. 6.3 From the fit, we find that A = 0.09, which 
corresponds to the attraction time, t = 0.05ms, of the system toward the critical 
solution, is similar across the three evolution variables. We further observe that the 
attraction time is similar to the departure time found in (47] which indicates that the 
eigenvalues of the decaying mode and the unstable mode have degenerate real parts. 
We note that these time scales are of an order of magnitude smaller than the cooling 
time scales of neutron stars reported in the astronomy and astrophysics literature 
[68], (69],|70j. Therefore, there is a high possibility that real astrophysical systems of 
cooling neutron stars pass through the threshold of critical collapse. Further, systems 
that undergo slow accretion and angular momentum loss with time scales larger than 
that reported here for the neutron star critical solution may also pass through the 
threshold and exhibit critical behavior. 

In order to understand the time scales of the oscillation phase of this critical so- 
lution, we now seek to compare the frequencies to that obtained via a perturbative 
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approach. We observe that the rest mass of the object undergoing critical collapse 
is higher than the maximum allowed mass for a TOV configuration with the same 
equation of state. Therefore, we consider a static non-rotating equilibrium TOV con- 
figuration having the same rest mass as the rest mass of the two TOV configurations 
producing the critical solution added together, namely at rest mass Mf, = 1.5. We 
recall from Chapter 5 that the critical set is a stationary or self-similar spacetime. 
In the case of axisymmetric neutron star systems considered in |47j . black holes form 
with a finite mass, and scale-invariance is not seen, therefore the critical collapse is 
of Type I, and the critical set is thus a stationary spacetime. In light of this, the 
critical set of this solution is non-radiative in principle, and we can then employ the 
tools presented in Chapter 4, ie. the 1=0 and even parity 1=1 mode perturbations, to 
determine the non-radiative pulsations of a corresponding non-rotating TOV config- 
uration. We obtain the finite-differenced TOV solution, namely the finite-differenced 
matter-energy density, the specific heat enthalpy, the pressure and the ADM mass 
function of the solution with respect to the isotropic radius of the configuration. 
We then manipulate these functions accordingly and insert them into Eq.s (4.70) to 
(4.72), and follow the rest of the procedure as mentioned in Section 4.3. Via this 
method, the fundamental I = frequency for a static equilibrium non-rotating TOV 
configuration is found to be u± p = 6.858 x 10 -3 . Similarly, for the even-parity I = 1 
mode, we insert the manipulated functions of the TOV solution into Eq. (4.75). The 
estimated fundamental frequency for this mode is found to be U2 P = 1.813 x 10~ 2 . 
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We note that uji p and U2 P are both one and two orders of magnitude smaller than 
uj\ and Co>2 respectively. This is a confirmation that the oscillations of the critical set 
of the solution are not caused by non-radiative pulsations of an equilibrium TOV 
configuration. There therefore exists the need of performing perturbative analysis on 
non-equilibrium rather than equilibrium background spacetimes in order to ascertain 
the nature of the oscillations exhibited by the neutron star critical solution. 

6.2 The neutron star critical gravitational collapse 
solution as a semi-attractor 

In [UJ , the critical index is seen to be similar across variations of different parameters 
of the initial data ie. the boost velocity, initial separation and equation of state. 
In line with the reasoning presented in Section 5.2, we acknowledge that via this 
method, universality cannot be strictly claimed for the neutron star critical solution. 
In this section, we present another method to prove universality. In this approach, we 
construct a neutron-star like initial data by solving the initial value problem where 
the matter field consists of two packets of matter whose densities are characterized by 
Gaussian distributions rotated about the axis of axisymmetry and follow a polytropic 
equation of state, p = np v e with k = 80 and T = 2. The Gaussian density distribution 
is as follows: 




(6.4) 
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Figure 6.4: The Gaussian matter density distributions are constructed so as to maintain 
the rest mass of the system as a constant. Therefore, when the Gaussian 
heights are increased, their widths are correspondingly decreased. Configura- 
tions 1 and 2 are less compact distributions whereas that of 3 and 4 are more 
compact. 
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Configuration 


M b 


Pc 


d 


v z 


NS 


1.6 


0.00056595 


27.6 


0.15 


1 


1.5 


0.00038698 


29.6 


0.1 


2 


1.5 


0.00039192 


29.6 


0.1 


3 


1.6 


0.00055501 


27.6 


0.12 


4 


1.6 


0.00064912 


27.6 


0.12 



Table 6.1: Mj, is the rest mass, p c is the neutron star central matter density /height of 
the Gaussian matter density distribution at t = 0, d is the center-to-center 
separation between the neutron stars/Gaussian packets, and v z is the boost 
velocity of the neutron stars/Gaussian packets along the z-direction of the 
grid. 

where A, B and C are constants that can be adjusted to yield different spacetimes. 
The normal coordinate condition (3 = and the geodesic slicing condition N = 1 are 
chosen for the initial data. The 3-metric 7 for the spacetime of this initial data is ap- 
proximately asymptotically flat. We shall comment further on how an approximately 
asymptotically flat spacetime differs from one that is exactly asymptotically flat in 
our analysis in Section 6.5. The tuning of the height and width of the Gaussian den- 
sity distributions produce a new 1-parameter family of initial data. The evolution of 
these initial data sets are carried out using the same numerical setup as used in [17] , 
namely, using the T-freezing shift and 1 + log slicing condition. The BSSN scheme is 
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also used in all our simulations as in |47j . 

We first set up this initial data set on the same numerical finite-differenced grid 
as that used in |47j . namely, on a thin slab of 323 points in the x and z-directions and 
5 points along the y direction. The Gaussian packets are then boosted to a head-on 
collision. The boost is performed via a coordinate transformation of the spacetime 
as follows: 



t' 




w 


v x w 


VyW 


v z w 




t 


x' 




v x w 




I (W-l)V X Vy ^ 


( (w-l)v x v z ^ 
\ v 2 ) 




X 


y' 




VyW 


i (w—l)v x v v \ 

\ V 2 > 




^ (w-l)v y v z j 




y 


zf 




v z w 


/ (w-l)v x v z \ 

\ V 2 > 


^(W-l)VyV Z ^ 


+ 




z 



(6.5) 



where v = (v x ,v y ,v z ) is the boost velocity vector, v 2 = v 2 , + Vy + v\ and w = 
(1 — f 2 ) -1 / 2 . Therefore, the values of the boost velocity vector components do not 
possess intrinsic physical meaning. However, the variation of the boost velocity is 
non-degenerate and well-behaved, ie. we can say that the increase or the decrease 
of the value of the boost velocity does indicate the isomorphic increase or decrease 
in the physical boost velocity of the objects in the system. The behavior of the 
axisymmetric object formed is observed. We note that this initial data is in a non- 
equilibrium state at the initial time, although it satisfies the constraint equations. 
The constraint equations with the two Gaussian distributions of matter and their 



boost velocities are solved using York's formulation (Section 3.4). Table ?? and 



Fig. |6.4| show the configurations that are constructed and their density distributions 
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along the direction of collision respectively. These initial data sets are distinct from 
the neutron star configuration in terms of their spacetime and matter properties, 
as they are clearly not obtained by varying different parameters of the neutron star 
configuration, nor are bound by the TOV spacetime and matter configuration. They 
are also different from each other in terms of their matter configurations. We examine 
the evolutions of these initial data sets to see whether they are attracted to the same 
oscillatory critical set of the neutron star solution. 

In order to clearly present this behavior of the evolution of these initial data 
sets, we choose several variables, namely the trace of the extrinsic curvature K at 
the coordinate (2.1,0.06,0.06), the z-component of the momentum density at the 
matter density contour of 7.5 x 10~ 3 , and the maximum matter density of the system. 
In Fig. 6.5, we draw the trajectories of these configurations using these variables. 
From this figure, we see that these variables correlate with each other such that they 
form bounded periodic orbits. In particular, we note that the correlation between 
the momentum density of a fluid element away from the center of collision and the 
matter density at the center of collision is analogous to that between the position and 
momentum of a simple pendulum. The oscillations of K before reaching the critical 
set is caused by the fact that the geodesic slicing is used at initial time. As the slicing 
switches from geodesic to 1 + log after the first time step, the value of the lapse at the 
grid boundary jumps discontinously from one value to another. We observe that these 
configurations are all attracted to the same critical set of the neutron star solution. 
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Figure 6.5: The trajectories of all the configurations start on the left and are attracted 
to the critical set on the right. The arrows point to the differing starting 
points of the trajectories. K denotes the trace of the extrinsic curvature 
at the coordinate (2.1,0.06,0.06) and S z denotes the z-component of the 3- 
momentum, Si, of the fluid element on 7.5 x 10~ 3 density contour, respectively. 
The natural logarithm is taken of the density at the center of the grid, namely 
Inpc, so as to greater distinguish the differing starting points of the trajectories 
in the phase space. 
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Figure 6.6: 

This provides evidence that the neutron star critical solution is a semi-attractor. 

6.3 Phase space analysis on the neutron star crit- 
ical gravitational collapse threshold 

In this section, we probe the structure of the critical solution in the framework of 
phase spaces. We note that Gaussian initial data sets possess additional degrees of 
freedom in the matter and in the spacetime as compared to the neutron star initial 
data. Using the additional degree of freedom in the matter configuration, we vary 
the boost velocity of the Gaussian packets together with the Gaussian heights and 
widths while maintaining the rest mass of the system at M b = 1.6389. In Fig. 6.6, we 
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Figure 6.7: Comparison between the neutron star and Gaussian packet critical surfaces 
in the rest mass-boost velocity phase space. 
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Figure 6.8: Comparison between the neutron star and the "corrected" Gaussian packet 
critical surfaces in the rest mass-boost velocity phase space. 
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Figure 6.9: 



plot boost velocity with respect to Gaussian height. In this phase space, we use the 
maximal slicing for the Gaussian packet initial data. We observe negligible deviation 
of the evolutions of these initial data sets with those of the Gaussian packet initial 
data sets in the previous section that were set with geodesic slicing. A main feature of 
the phase space depicted in Fig. 6.6 is that there is a turning point at approximately 
(p c ,f) = (0.0009,0.12). This is an indication that for a certain compactness, the 
Gaussian packet system crosses two phase thresholds as the boost velocity is increased, 
ie. from a neutron star-like phase to a black hole phase and back to a neutron star-like 
phase. The cross from the neutron star-like phase to the black hole phase is expected 
as the boost velocity increases entails the increase of gravitational energy. However, 
we note that the cross from the black hole phase back to the neutron star-like phase 
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when the boost velocity is increased further, is a surprising physical behavior. 

Given the limitation imposed by the grid size, we are unable to extend the curve 
further to the left of the phase space. An initial data configuration to the left of the 
critical surface collapses to a black hole and one to the right results in a neutron star- 
like object. The extent of this critical surface is indicative of the size of the attraction 
basin of the neutron star critical set with rest mass 1.6389. Beyond the right edge 
of the critical surface at this rest mass, no critical collapse behavior is observed, ie. 
configurations with Gaussian packets with compactness p c > 0.0009 are not attracted 
to the critical surface at all. To further explore the extent of the attraction basin, we 
construct another phase space using the boost velocity and the separation distance 



between the Gaussian packets. Fig. 6.9 shows two critical surfaces of which the left is 
where the system at a fixed separation distance passes the threshold from the neutron 
star-like phase into the black hole phase and passes back into the neutron star-like 
phase through the right critical surface. Both the critical surfaces end at d/2 ~ 14 
at the top of the phase space. The bottom extent of the right critical surface is 
limited when the separation distance between the two Gaussian packets decreases 
to the point that the packets merge to become one single packet. As the total rest 
mass of this single packet, ie. 1.6389, is more than the maximum rest mass of an 
equilibrium neutron star configuration, ie. 1.6087 for configurations with T = 2.0, 
the system collapses into a black hole even when there is no implosion velocity. We 
thus set the rest mass of the system to be 1.6, which is below the maximum rest mass. 
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As the maximum point of the rest mass-central density relation for equilibrium star 
configurations (Fig. 4.1) is not attracted to the critical set, we study the transition 
of the neutron star critical set attraction basin by increasing the central density of 
the Gaussian packets so that it approaches the maximum point, ie. p c ~ 4.0 x 10 -3 . 
We find that the oscillations of the critical set decrease in amplitude. The average 
of the oscillations in the central lapse function N shifts up whilst that of the central 
density shifts down. At p c ~ 1.833 x 10~ 3 , the threshold between the black hole 
and neutron star-like phases disappear. At central densities beyond this value, the 
system oscillates at the normal mode frequencies of its corresponding equilibrium 
star configuration not only when there is no implosion velocity but at all other values 
of the implosion velocity up till 0.8. However, we observe that the transition point 
across p c ~ 1.833 x 10~ 3 itself is characterized by a critical set between the black 
hole phase and the neutron star-like phase, namely, when the central density is more 
than this threshold value, the single Gaussian packet settles into an oscillating state 
about a corresponding equilibrium star configuration, while it collapses into a black 
hole when the central density is less. 

We also construct a phase space of the boost velocity and the rest mass to compare 



with the critical surface found with that for the neutron star initial data (Fig. 6.7). 
We note that there is also a turning point behavior in both critical surfaces, similar 
to that found in the boost velocity-Gaussian height phase space above. However, 



another striking feature of Fig. 6/7 is that the Gaussian packet critical surface is 
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shifted to the bottom of the neutron star critical surface due to the fact that the 
spacetime for the Gaussian packet initial data is only approximately asymptotically 
flat, ie. the asymptotic value of the 3-metric is slightly less than 1. This difference 
in the spacetime causes the rest mass summed up over each hypervolume element of 
the Gaussian packet initial data, given by: 



M b = j ^pWd 6 x, (6.6) 

to be less by a scale with respect to that for the neutron star initial data. However, 
based on the similarity of their evolution trajectories, we emphasize that the critical 
solution represented by the Gaussian packet initial data point with Mb = 1.5505 is 
the same as that represented by the neutron star initial data point with Mj, = 1.6379, 
Gaussian Mb = 1.5600 the same as neutron star Mb = 1.6382, and so on. This 
represents a scaling effect in the critical surface of the rest mass-boost velocity phase 
space caused by the difference in spacetime as mentioned earlier in this paragraph. 
We also vary the spacetime by changing the height and width of the metric function 



in the region occupied by the Gaussian packets. The black boxes in Fig. 6.7 represent 
critical solutions found when the metric function is varied in such a manner. 

To test universality, we set the spacetime of the Gaussian packet initial data 
back to asymptotic flatness. The Gaussian critical surface in this phase space shifts 
back up approaching that of the neutron star as the asymptotic behavior of the 
Gaussian packet initial data spacetime approaches that of the neutron star initial 
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data (Fig. 6.8). We see a common rest mass range but not a common boost velocity 
range. The discrepancy in the boost velocity ranges is due to the fact that the 
Gaussian packet initial data possess different coordinate systems than the neutron 
star initial data. We observe that a critical surface with a common boost velocity 
range can also be obtained, indicating that the coordinate system of the initial data 
can be freely adjusted using the additional degrees of freedom in the Gaussian packet 
configuration. In fact, we are able to obtain a family of critical surfaces in the rest 
mass-boost velocity phase space with a common rest mass range and within the boost 



velocity range indicated in Fig. 6.6 for the rest mass of 1.6389, ie. between 0.06 and 
0.21. As rest mass is conserved, the overlap in the rest mass range of both critical 
surfaces confirms that neutron star-like systems evolve toward the same critical sets 
as the neutron star system. 

Using the ADM mass in replacement of the boost velocity, we find similar critical 



surfaces with turning points (Fig. 6.10 and Fig. 6.11). The ADM mass is calculated 



in our numerical simulations using a volume integral as follows: 

M ADM = [ ¥>[ Pe + ^riAijAv _ 2 K 2 _ ^-4)] y~ rf 3 x? (6 7) 

J 167T 6 

using the conformal variables that have been defined in Section 3.2 and the matter- 
energy density as mentioned in Section 4.1. In line with the definition that is also 
given in Section 4.1, this mass is the total gravitational energy of a system that 
is isolated at spatial infinity, and requires a quasi-isotropic gauge choice in which 
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jijj ~ Oir 3 ) as r — > oo [85J, a condition that is instantaneously realized by the 
initial data before it gravitationally radiates, even though in non-equlibrium states. 



In Fig. 6.10 we also compare the rest mass-ADM mass critical surface between 
that of the Gaussian system and that of the neutron star system. Branch 1 and 
branch 2 of the critical surface in the rest mass-ADM mass phase space represent two 
different phase thresholds. In Section 6.5, we comment on the indices of these critical 
solutions and their comparisons between the Gaussian system and the neutron star 
system. 

Results of this section thus indicate that the neutron star critical solution is uni- 
versal with respect to spacetime and matter configuration and reinforces the claim 
in the previous section that the neutron star critical solution is a semi-attractor. We 
will further support this evidence by analyzing the properties of their critical indices 
in Section 6.5. 

6.4 Spacetime and matter properties of the critical 
solution 

In this section, we explore the various spacetime and matter properties of the neu- 
tron star critical solution. For this purpose, we first ensure that our critical collapse 
simulations are converging with respect to the finite-differencing resolution up till 



the duration of the critical solution (refer to Appendix A). Fig. 6.12 shows the first- 
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Figure 6.12: 



order convergence in the Hamiltonian constraint violation of our simulations using 
the GRAstro-2D code at t — 243M , which is well into the fourth oscillation of the 
matter density and 4-scalar curvature of the critical solution at the center of colli- 
sion. The anomaly at z = is an expected effect caused by the truncation in the 
finite-differencing at the center of symmetry which has been shifted (to z = 0.03 for 
dx = 0.06, z = 0.06 for dx = 0.12, and so on) by a staggered grid. We then measure 
the baryonic mass enclosed within different density contours throughout the merged 
object, the proper polar circumferences of these density contours, and their equatorial 
and proper radii, and observe their evolution throughout the critical set oscillation 
phase. The measure of the baryonic mass enclosed between different density contours 
throughout the merged object has the advantage of being a fully-geometric measure. 



112 



Chapter 6 



Critical gravitational collapse of neutron stars 



0.45 
0.40 
0.35 
0.30 
0.25 
0.20 
0.15 
0.10 
0.05 
0.00 

80 



$5 70 

e 

C 60 
| 50 

^ 30 




H 1 1 1 1 h 




H 1 1 1 1 h 




2.5e-4 
5.e-4 
7.5e-4 
2.5e-3 
3.5e-3 
5.e-3 
7.5e-3 



120 140 160 180 200 220 

coordinate time, t 



Figure 6.13: Right legend shows contour densities at which measures are performed. Polar 
and equatorial proper radii shown on the right are measured at the 5.0 x 1CT 4 
density contour. The polar direction is taken to be the x-axis on the xz-plane 
of the grid and the equatorial direction is taken to be the z-axis on the same 
plane. The equatorial direction is thus along the direction of collision of the 
neutron stars. 
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In Fig. 6.13, we observe that this measure exhibits oscillations that are correlated be- 
tween different regions of the merged object, whereas the corresponding proper polar 
circumferences do not exhibit such well-defined oscillations throughout the critical set 
dynamics, ie. on the time scale of t ~ 40. The well-defined oscillations in the bary- 
onic masses summed up between the density contours thus indicate the existence of 
matter fluxes or a slushing back and forth of the matter through the different density 
contours inside the object and are not due to any correlated changes in contour cir- 
cumferences. The proper radii however exhibit prolonged oscillations only within the 



region bounded by a specific density contour. Fig. 6.14 shows the oscillations of these 
radii outside this region exhibiting damping into a stationary state. These oscillations 
are correlated between the polar and equatorial directions ie. the polar proper radius 
is in phase with the oscillation of the central density while the equatorial radius is 
out of phase. However, we note that the proper radii oscillations presented here have 
not eliminated the effect of the lapse at the outer limit of the proper radii integration. 
A good method to isolate this effect is to set a diagnostic spacetime with the normal 
coordinate choice and where the 11-component of the 3- metric is set to 1. With such a 
spacetime, the hypersurface will always be perpendicular to the timelike normal unit 
vector n and the proper radius oscillation will reflect the real oscillation of the den- 
sity contour of the neutron star [SH] ■ In order to gauge the stationarity of the merged 
object, we also measure the power of gravitational radiation emission throughout the 



dynamical process leading to and that of the critical set itself. Fig. 6.15 shows the 
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power that is calculated using the quadrupole formula ([12], [25]), namely: 

dF " 5 Q ^' Q ^ ' (6 - 8) 

where Qy = J p{x % x^ — ^5 l:i r 2 )d 3 x. From this figure we see that the energy radiated 
gradually decreases as the system approaches the critical set at t ~ 110, which is the 
approximate time of merge of the two neutron stars. The power remains at values very 
close to zero until t ~ 165 when it jumps up to positive values that increase with time. 
We note that these sharp jumps may be due to noise created by finite-differencing 
the time derivatives of Qij |55j. The average of the power throughout the critical set 
dynamics remains very close to zero. This suggests that the oscillations of the critical 
set may not damp out with time thus pointing to the possibility that the critical set 
is a limit cycle rather than a limit point. Also, it shows that at t ~ 110, when the 
oscillation phase of the critical solution starts, the system has already settled down 
to a configuration that is stationary enough and therefore very near to an equilibrium 
state. 



6.5 Properties of the critical index 

As mentioned in Section 6.3, this section will contain further elaboration of the prop- 
erties of the critical indices of the different Gaussian packet initial data configurations 
as compared to that of the neutron star. 

Analogous to Eq. (5.6), the critical index is calculated using two departure thresh- 
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olds 5%, 20%, where the departure threshold is taken as the time when the percent 
difference between the evolution variable of a solution leaving the critical set and that 
of the critical solution reaches a certain value [UJ . In this case, we take the evolution 
of the lapse at the center of collision. We denote the departure thresholds as T 05 , T 2 
and so on according to the percent difference used. For each critical point, we find 
that T p is linearly proportionate to In \p — p * |, where p is any initial data parameter 



that we choose to vary while fixing all others. In the case of Fig. 6.16, the initial 
data parameter varied is the boost velocity of the Gaussian packets. Straight best fit 
lines are therefore obtained and the standard deviations of their slopes are plotted 
as errorbars in this figure. We note that this errorbar only represents one type of 
error, ie. the error affecting how the evolution trajectories follow the power-law be- 
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havior characteristic of critical and near-critical solutions. Other errors that affect the 
numerical results include the error caused by the finite differencing scheme and the 
error affecting how the evolution trajectories leave the critical set in an exponential 
manner. 

We check the convergence of the critical index with respect to the finite-differencing 



resolution dx. Fig. 6.17 shows the critical indices for different departure thresholds at 



three different resolutions for a sample Gaussian system. At the departure threshold 



of 0.2, Fig. |6.18| shows an approximate first order convergence with respect to the 
resolution. The grid boundary used for the simulations with resolution dx = 0.16, 
namely, 41.0, is slightly larger than that with dx = 0.12 and dx = 0.24, namely 
38.5, due to the restrictions on the number of grid points imposed by the multigrid 
initial value problem solver, which accounts for the slightly larger than first order 
convergence demonstrated by the dx = 0.16 point in this figure. 



Fig. 6.19 shows the rest mass dependence of the Gaussian packet critical index 
compared with that of the neutron star, where we observe that the critical index for 
both the Gaussian packet and neutron star initial data increases with rest mass of the 
system. This indicates that the time scale of departure of a near-critical solution from 
the critical set decreases as the rest mass of the system increases. The overlapping of 
the critical indices of both the Gaussian and the neutron star systems shows that the 
neutron star critical set is universal with respect to different initial data, where for 
the same rest mass and within a small interval of ADM masses, the Gaussian initial 
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data produces critical indices that are within the errorbars of those produced by the 
neutron star initial data. The universality of the neutron star critical index with 
respect to these different initial data configurations is a crucial point in establishing 
the quality of the neutron star critical set as a semi-attractor. 

To answer the question of the parameter dependence of the critical index, we 
consider the critical solutions represented by each point in the ADM mass-Gaussian 
height phase space and study how the critical indices change along this critical surface. 



Fig. |6.20| shows the variation of critical index along this critical surface with respect 
to the ADM mass fixing the rest mass at 1.6389. The solid line is a fit of the average 
of the indices obtained at departure thresholds 0.05 and 0.2. We note that the first 
order convergence of these critical indices with respect to resolution holds for different 
departure threshold regimes for configurations with different ADM masses. This is 
due to the fact that for configurations with higher compactness, and hence lower 
ADM mass, the grid boundary is further away from the physical system than for 
configurations with lower compactness, or higher ADM mass. Convergence is better 
achieved when the grid boundary is further away from the physical system as there 
is less interference in the spacetime due to gravitational radiation. 



In Fig. 6.10 we take the rest mass of 1.65 and extract the critical index for 



the critical points on both branches 1 and 2 of the critical surface. In accordance 



with Fig. 6.20 we also find a discrepancy between these critical indices that goes 



beyond their errorbars, indicating that critical points on branches 1 and 2 of the 

121 



Chapter 6 



Critical gravitational collapse of neutron stars 



critical surface represent two different phase thresholds respectively. On branch 1 
of the critical surface with this rest mass, the critical index obtained is 10.9 ± 0.1 
at departure threshold 0.05 and 10.84 ± 0.08 at departure threshold 0.2, whilst on 
branch 2, the critical index is found to be 10.26 ± 0.06 at departure threshold 0.05 
and 10.24 ± 0.07 at departure threshold 0.2. 

As mentioned in Section 6.3, the ADM mass of the initial data is evaluated using 



Eq. (6.7). The constancy of the ADM mass guarantees that Einstein's field equa- 
tions for the system are satisfied for all time [33]. Therefore, it is by definition a 
conserved quantity. However, when evaluated at a finite boundary for a system that 
gravitationally radiates, it is seen to decrease with time as gravitational energy is 
carried out of the finite computational grid by gravitational radiation. The mass that 



is evaluated using Eq. (6.7) in such a scenario does not have strict physical mean- 
ing although it can be used in an approximate conservation equation primarily to 
monitor the accuracy of 3 + 1 numerical simulations |73j . The fact that the quantity 
lacks a strict physical meaning is because the hypersurface is not spatially isolated 
at a time after gravitational radiation. This is evident in that no matter how far 
one pushes the computational boundary spatially outward, one would encounter the 
gravitational radiation that is propagating to null infinity. In other words, the ADM 
mass measurement does not by construction separate the gravitational energy of the 
physical system that we are interested in and that of gravitational radiation. Even 



for Eq. (6.7) to hold approximately in the computational grid, the boundary of the 
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computational grid has to be asymptotically flat and thus free from the presence of 
any gravitational radiation, the ambient spacetime has to be fully stationary and 
fulfil the right gauge conditions as mentioned earlier. Therefore, in order to employ 
this equation even in an approximate sense to evaluate the gravitational energy of 
a physical system after gravitational radiation leaves it into the ambient spacetime, 
one has to do the measurement at a very late time of the critical set dynamics when 
one is sure that gravitational radiation is no longer present inside the computational 
grid and is no longer interacting with the computational grid boundary in any way. 
However, none of the initial data in numerical simulations can be set directly on the 
critical threshold in principle, implying that they will either evolve into a neutron 
star-like object or collapse into a black hole at one point in time. Also, due to the 
accumulation of numerical errors, in practice, we are not able to accurately simulate 
the critical set dynamics until very late time. Intrinsically as well, even though by 
definition stationary and an initial data that dwells directly on the threshold will 
inevitably evolve into it, the critical set itself is not a TOV equilibrium configura- 
tion. A measurement that is based on a strict physical meaning however is the Bondi 
mass measurement which evaluates the gravitational mass of the system isolated at 
null infinity. This measurement puts any gravitational radiation to the past of the 
null surface of the evolving physical system as the radiation itself propagates to null 
infinity. 

Considering the question of the parameter dependence of the critical index on the 
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above-mentioned conserved geometric quantities, we are basically asking a question of 
how many classes of neutron star semi-attractors exist according to the rest mass and 
ADM mass. A natural approach to this question that we have employed is to vary the 
rest mass and when the rest mass is fixed, to vary the ADM mass of the initial data 
or the intrinsic gravitational energy of the physical system minus the energy that will 
be radiated away during its evolution to the critical set. We recall that the neutron 
star-like initial data that we have constructed possess additional degrees of freedom 
in comparison to the neutron star initial data, that enables us to fix the rest mass and 
vary their intrinsic gravitational energies. This intrinsic gravitational energy consists 
of the binding energy inherent in the Gaussian packet system and in principle will not 
be radiated away unless there is a physical mechanism that diffuses this energy, eg. 
neutrino emissions, which is absent in our numerical simulations. Therefore, when 
we vary the compactness of the Gaussian packets, we are intrinsically providing dif- 
ferent gravitational energies to these systems which will be conserved in time. This 
is evidenced in the fact that when we fix the boost velocity of the initial data sets 
and vary only the height of the Gaussian matter distribution, and correspondingly 
the compactness of the matter packets, by Ap c ~ 2 x 10~ 4 , we see a change in the 
ADM mass by 0.01. This change is free from the factor of the mechanical energy in- 
herent in the boost velocity that can be radiated away by the gravitational radiation 
- a physical system does indeed radiate more when it is moving faster. If they are 
radiated away by some physical mechanism which we can choose to impose on our 
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numerical simulations, or by numerical viscosity, and if we assume that the gravita- 
tional energies of the remnants which in our case is the semi-attractors, are all the 
same, our entire investigation of whether there exists a 2-parameter dependence of 
the critical index becomes a meaningless endeavor, as the motivation of the question 
depends on the premise that we are indeed able to manipulate the additional degrees 
of freedom in our initial data sets to vary the compactness and correspondingly the 
intrinsic gravitational energies of the physical system. In other words, the premise 
of our question in the first place is that the semi-attractors themselves can have dif- 
fering gravitational energies for the same rest mass. Also, the critical index itself 
is a geometric quantity and is evaluated during the semi-attractor dynamics. The 
observation that they indeed differ even when we adopt the above assumption that 
the gravitational energies of the semi-attractors are all the same, would result in the 
observation that semi-attractors with the same gravitational energies have differing 
critical indices, a confounding observation that not only is unlikely, but throws this 
entire question of 2-parameter dependence meaningless. Considering realistic astro- 
physical systems where dissipative physical mechanisms which depletes the binding 
energies of the packets of matter do indeed occur, it is still very unlikely that a fine- 
tuned balance exists for the gravitational dynamics, where the gravitational energies 
are all depleted by nonlinear physical mechanisms in just such a way that they all end 
up the same for all the remnants. We cannot rule out such a scenario, which would be 
an interesting physical property for neutron star semi-attractors, but such a scenario 
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Figure 6.21: 



would still violate the premise of the 2-parameter dependence question as mentioned 
above. In addition, we have to recall that the neutron star semi-attractor is not a 
TOV equilibrium configuration, which does not allow us to draw strict analogies with 
stable equilibrium TOV configurations which has a unique ADM mass for each of its 
rest masses. 



Fig. 6.21 shows the dependence of the critical indices with respect to the averages 



of the gravitational masses, (Mg), measured by Eq. (6.7) on the finite computational 
grid over the period of two oscillations in the lapse function at the center of collision 



of each configuration along the critical surface in Fig. 6.11 Fig. 6.22 shows the evolu- 



tion of this gravitational mass calculated using Eq. (6.7) on the finite computational 



grid, which we denote as just Mq, for a sample Gaussian packet system that evolves 
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toward the critical set and dwells there until t ~ 240 for two differing resolutions and 
differing grid sizes. The observation that the time shift in the curves for the simula- 
tion with the boundary at x = z = 41.0 from the one at x = z = 38.5, namely t ~ 3 
measured at the first dips of the curves, which is comparable to the spatial shift using 
c = 1, and that the oscillations do not converge away with resolution, shows that the 
oscillations in M G is caused by the interaction of gravitational radiation with the com- 
putational grid boundary. As mentioned above, this renders the measurement itself 
not only numerically inaccurate but without even an approximate physical meaning. 
Nonetheless, we still see a dependence of the critical index with respect to (M G ). As 
mentioned above, the fact that the critical indices do vary when we vary the com- 
pactness of the Gaussian packet systems, will not change whether we plot them with 
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respect to either the ADM masses or with (Mq), unless all the (Ma) values - assuming 
that they are even accurate and possess an approximate physical meaning - end up 
being the same up to numerical error, which is both physically unlikely and a viola- 
tion of the premise of the question that we sought to answer. Therefore, even though 
we do not rule out the possibility that the premise of our question is unfounded, ie. 
that indeed, neutron star semi-attractors cannot in principle, even though they are 
not in a TOV equilibrium state, intrinsically possess differing gravitational masses 



for the same rest mass, we conjecture from our observations (Fig. 6.20 and Fig. 6.21) 
that there is still a high likelihood that there exist classes of semi-attractors that can 
be labeled by both the rest mass and the gravitational mass of the neutron star-like 
system, due to the fact it would be a very rare although interesting occurence that 
there exists a fine-tuned balance in the gravitational dynamics such that gravitational 
energies of physical systems with intrinsically different binding energies are dissipated 
by nonlinear physical mechanisms in the matter to end up at the same value. In addi- 
tion, that there exist such classes of semi-attractors is also in line with our observation 



of two distinct phase thresholds for the critical surface in Fig. |6.10| at the same rest 
mass, where the distinct phase thresholds represent two differing semi-attractors with 
the same rest mass but with differing ADM masses. 



Fig. 6.23 shows the convergence of the critical indices for two configurations on 



the ADM mass-Gaussian height phase space, ie. the configurations with Madm 



1.553 and Madm = 1-5868 respectively, whereas Fig. |6.24| shows the convergence 
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of the rest mass averaged over the period of the first two oscillations of the lapse 



function at the center of collision, which we will denote as (Mb), and Fig. 6.25 the 
convergence of (Mg)- We see that (Mb) and (Mg) converge on a higher order for the 
configuration with higher compactness and vice versa, ie. for the configuration with 
higher compactness, the (Mb) and (Mg) converge in a fourth order manner and fifth 
order manner respectively, but for the less compact configuration, they both converge 
in a second order manner. However, the (Mb) for both configurations converge to 
1.64 up to the numerical error caused by the finite differencing, whilst their (Mg) 
values converge to different values on the order of 0.01 respectively, ie. 1.38 for the 
more compact configuration and 1.41 for the less compact configuration, outside the 
error bars from the finite differencing. The convergence of the (Mb) to the same 
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value is expected due to the fact that the rest mass of the initial data have been set 
to the same up to the eleventh decimal, and the fact that rest mass is a conserved 
quantity even when evaluated on a finite computational grid. The variations in the 
twelfth and above decimals of the rest masses of the different initial data or in their 
{Mb) values in the limit of infinite resolution are completely random. As for the 
critical indices, similar to (M G ), their values for both configurations are observed 
to converge to different values, where the more compact configuration exhibits first 
order convergence at a faster rate than the one for the less compact configuration, 
which is to be expected due to the different distances of the computational grid 
boundaries from the physical systems as mentioned previously. These observations 
imply that in the limit of infinite resolution, where a large portion of numerical 
errors is eliminated, even though putting aside the ambiguous physical meaning of 
(Mq), the critical index depends on both the compactness, with its corresponding 
binding energy present in the physical system, as well as the rest mass of the system, 
which is in line with the conjecture we make that is mentioned earlier in this section. 
Since we have seen that convergence of the critical index is achieved at an even 
faster rate when the computational grid boundary is located further away from the 
physical system, when the boundary of the computational grid is pushed further out 
for each configuration, we expect to obtain similar observations that point to the 
above-mentioned conjecture. 
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7.1 Conclusions 

In this chapter, we review the results obtained in this dissertation in accordance with 
the objectives we had in mind. In order to understand the time scales involved in non- 
rotating neutron star critical collapses, we first study the time scales of non-radiative 
pulsation modes of single static neutron stars. We find that these time scales are an 
order of magnitude higher than the time scales involved in the non-rotating neutron 
star critical collapses observed in simulations. We thus confirm that the pulsations 
exhibited by the quasi-stationary object during the critical collapse dynamics are not 
that of non-radiative pulsation modes of single static neutron stars. This points to 
the fact that the quasi-stationary object is not in a state of TOV equilibrium. 

Via curve fitting and Fourier transform, we determine the time scale of attraction 
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of the neutron star system toward the neutron star critical set. This time scale is 
an order of magnitude smaller than that reported of neutron star systems under- 
going slow cooling. We thus conjecture that dynamical neutron star systems with 
significant amount of kinetic energy undergoing slow cooling via formation processes 
and radiation transport or systems undergoing slow accretion and a gradual loss of 
angular momentum have a possibility of undergoing critical collapse. 

By constructing new neutron star-like initial data not restricted by the TOV 
equations and studying the behavior of their evolutions, we confirm that the neutron 
star critical set is a semi-attractor. We employ the freedom in the gauge choice and 
the matter configuration to construct different such neutron star-like initial data and 
by measuring their spacetime and matter 'flows', we observe that these initial data 
evolve towards the neutron star critical set. The critical indices extracted from these 
evolutions are the same up to numerical error as that obtained from neutron star 
simulations with the same rest mass and ADM mass. We thus reinforce the claim 
that the neutron star critical set is universal with respect to initial data, going beyond 
the perturbative regime of the critical set. 

We further study the properties of the neutron star critical set dynamics itself. We 
observe that the oscillations of matter contained within different density thresholds 
of the object undergoing critical collapse exhibit phase differences across a certain 
density threshold, indicating the existence of matter fluxes across a certain density 
threshold within the object. We also find that there exists an outer envelope within 
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the object that damps out at a rate higher than the central region. In addition, 
the proper circumferences of different density thresholds within object are nearly- 
constant. To further shed light on whether the semi-attractor is a limit point or a 
limit cycle, we measure the power of the gravitational radiation emission throughout 
the critical collapse dynamics. It is observed that the average of this power during 
the critical collapse process is very close to zero, suggesting that the oscillations of 
the critical collapse object may not damp out with time and thus that the semi- 
attractor is a limit cycle. Combined with the possibility of realistic neutron star 
systems undergoing critical collapse based on the time scale study mentioned above, 
this characteristic of the neutron star critical solution, which is characterized mainly 
by the existence of a universal unstable mode and its specific gravitational radiation 
signature, may cause a queiscent effect on the gravitational radiation spectrum of 
realistic neutron star systems. 

Using the additional degrees of freedom available to us in the new neutron star-like 
initial data we construct, we also investigate the properties of the characteristic time 
scale of the neutron star critical set, which is the time scale of its unstable mode, 
describing the time scale of departure of a near-critical evolution from the critical 
set. We find that this critical index depends almost linearly on the rest mass of the 
system. We also observe that there is a possibility of a 2-parameter dependence of 
the critical index, where the parameters can be taken as the rest mass and the ADM 
mass of the system, both of which are conserved quantities. This indicates that there 
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exist classes of neutron star semi-attractors that can be labeled by both the rest 
mass and the gravitational mass of the system. This finding extends the observation 
that for a given rest mass, there exist two phase thresholds through one of which the 
system passes from the neutron star to the black hole phase at a lower ADM mass 
and through the other of which the system passes from the black hole phase back 
to the neutron phase at a higher ADM mass. We note that the observation of the 
existence of two phase thresholds for a given rest mass is similar to that obtained 
from the neutron star initial data. 

We then analyse the phase space properties of the neutron star critical set, by 
constructing phase spaces using the parameters of the rest mass, ADM mass and 
central density of the new system. For both the rest mass-ADM mass and ADM mass- 
central density phase spaces, we find turning point features with the aforementioned 
two-threshold feature. We note that the new system provides the additional degrees 
of freedom to construct the ADM mass-central density phase space. The extent of the 
critical surfaces in these phase spaces gives an indication of the size of the attraction 
basin of the neutron star critical set. For a certain rest mass, we determined the 
range of central densities and ADM masses of configurations that exhibit critical 
gravitational collapses, which helps us draw an approximation of the range of realistic 
neutron star or neutron star-like configurations that can exhibit critical phenomena. 

We further explore the boundaries of the attraction basin of the neutron star 
critical set by constructing a phase space using the separation distance and the boost 
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velocity (which also corresponds to the ADM mass) of the new system. In decreasing 
the separation distance, we move toward an initial configuration that consists of a 
single neutron star-like object with a varying implosion velocity Fixing the rest 
mass at a slightly lower value than the maximum rest mass allowable for a single 
equilibrium TOV configuration with T = 2, we find that there exists a boundary 
whereby the threshold disappears when we increase the central density towards the 
maximum rest mass point. 

7.2 Future work 

In order to further clarify the implications of the results of this dissertation on realistic 
astrophysical observations, we will extend the analysis on neutron star systems that 
possess angular momentum. To do this we intend to construct a similar neutron star- 
like initial data but with the additional incorporation of angular momentum, in both 
the irrotational and corotational orientations in accordance with the spin orientations 
that have so far been observed in realistic neutron star systems. This will also help 
us to see if an object that is simultaneously rotating and oscillating can undergo a 
dynamical process whereby the critical set is characterized by a limit cycle, which 
poses an interesting scenario in the theory of general relativity itself. 
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A The finite-differencing scheme 

The finite-differencing scheme is a numerical scheme where a grid is introduced as 
the domain of functions. Values of functions are evaluated on points on the grid 
thereby discretizing the functions according to the resolution of the grid. Using 
Taylor's expansion, an unknown function u(x, y) on a 2-dimensional grid is evaluated 
as follows: 

du ( Ax\^ d^u 
u(x + Ax, y ) = u(x , y ) + Ax—(x, y ) H — (x, y ) + 

where: 



du ( Ax^ d^u 

u{x + Ax, y ) - u{x - Ax, y ) = 2Ax— + + 0[(Ax) 5 }. (A-2) 



The partial derivative of u(x,y) with respect to x can thus be evaluated using a 
second-order central operator as: 

Ou _ U i+ ij - Ui-ij u 2 



dx 2 Ax 
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where i,j are the indices of the grid points in the x,y directions respectively. In turn, 
the next higher derivative can be written as: 

d 2 u _ u i+1 - 2u i:j + Ui-tj 



dx 2 (Ax) 



+ 0[(Ax) 2 ]. (A-4) 



Eq.s (A-3) and (A-4) describe a second-order finite-differencing scheme, as the dom- 
inant term on the right is second-order with respect to the grid resolution Ax. In 
addition, the order of the dominant term is preserved as we go to higher orders of 
the derivative. The scheme is thus said to be nth-order following the order of the 
dominant term of partial derivatives with respect to the grid resolution. Convergence 
of finite-differencing schemes is reached when Ax — > 0. At this limit, functions be- 
come in principle infinitely accurate. Without prior knowledge of a function, we are 
still able to check how accurately it is evaluated by checking its convergence with 
respect to the grid resolution. In numerical general relativistic hydrodynamics codes 
implementing finite-difference schemes based on the 3+1 formalism, convergence is 
checked by performing a simulation using different grid resolutions and observing how 
the violations of the Hamiltonian Eq. (2.37) and momentum constraint Eq. (2.38) 
scale with respect to the grid resolution. 



B Push-forward of vectors and pull-back of forms 

The push-forward of a vector along another vector field congruence is defined as 
an operator that acts on the vector such that it produces a new vector whose Lie 
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derivative along the vector field congruence is zero. Consider a vector B being pushed 
forward along the vector field congruence A. Let t be the affine parameter along A. 
The Lie derivative of B along A is as follows: 

£ A B = [ A, B] = Jim B ~ , (B-5) 

where ^ t B is the push-forward of B along A. The push-forward operator thus trans- 
ports the vector B along the vector congruence A such that the new vector *0 t B forms 
a coordinate basis with the form dual to A. Given this definition, we see that the 
push-forward of a vector acting on a form is the same as the vector acting on the 
pull-back of the form. The push-forward of a vector can also be seen as an active 
transformation on the vector whereas the pull-back of a form can be seen as a pas- 
sive transformation on the form while actively transforming the coordinate system in 
which the form is defined. 

C Geometric units 

The gravitational constant G has the units of [L 3 T~ 2 M~ 1 ], where L denotes the 
dimension of length, T denotes the dimension of time and M denotes the dimension 
of mass. Using the mass of the sun in metric units as a standard and the value of the 
constant as well as the speed of light in the same metric units, the unit for G can be 
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set dimensionless as follows: 

L 3 T~ 2 M~ l = 6.67 x KT 11 
(L/TfL = 1.33 x 10 20 

L = 1.33 x 10 20 /(3 x 10 8 ) 2 

L = 1.476 x 10 3 m, (C-6) 

which denotes 1M of length as equivalent to 1.476 x 10 3 m. From this relation, we 
can similarly write all other variables in terms of M , such as: 

lM Q oftime = 4.92 x 10~ 6 s (C-7) 

lM Q ofmass = 1.989 x 10 30 kg (C-8) 

lM Q of energy = 1.786 x 10 67 J. (C-9) 
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